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A spectral  estimator,  called  the  Wigner  distribution 
(WD),  is  investigated.  The  WD  in  the  context  of  the 
time-frequency  representations  of  signals  and  its  relation 
to  alternative  t-f  representations  is  examined.  The 
relationship  between  the  WD  as  a simultaneous  t-f  signal 
representation  and  the  one-dimensional  time  and  frequency 
domains  of  deterministic  signals  forms  an  important  part  of 
this  discussion.  The  effects  of  moving  averages  performed 
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on  the  WD  are  considered  in  terms  of  windowing  and 
filtering.  Window  and  filter  selections  for  maximum  signal 
concentration  are  determined. 

The  WD  is  applicable  to  analysis  on  continuous 
signals.  The  discrete-WD  (DWD)  is  the  discrete  analog. 
The  properties  of  the  DWD  and  its  relation  to  a digital 
signal  processing  (DSP)  setting  are  investigated.  A 
detailed  comparison  between  the  DWD  and  the  conventional 
FFT  based  moving  window  method  is  examined  both 
experimentally  and  theoretically  on  signal  types  whose 
dispersion  index  is  large.  The  smoothed-pseudo-DWD  is 
related  to  filtering  and  modulation  operations  on  the 
signal  then  experimentally  shown  to  improve  estimation  for 
the  two  primary  conditions  under  which  the  DWD  produces 
poor  results:  multi-component  signals  and  low 

signal-to-noise  ratios.  These  results  are  applied  to 
estimation  of  a multi-path  radio  example  in  the  presence  of 
severe  receiver  noise.  Emphasis  is  placed  on  the  DWD  as  a 
DSP  tool  and  the  theoretical  development  of  the  DWD  is 
heavily  supplemented  with  physical  examples  and 
applications  directly  related  to  signal  processing  and 
analysis.  Where  appropriate,  physical  discussions  replace 
or  supplement  formal  mathematical  descriptions. 

Wigner  processor  configurations  tailored  to  signal 
type  and  system  application  are  provided.  General  Wigner 
processing  using  standard  FFTs  are  addressed  and  a means  of 
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computing  two  N point  DWD  time  slices  via  one  N/2  point  FFT 
is  applied  to  the  case  of  real  signals.  A high-speed  DWD 
(and  spDWD)  processor  for  complex  data  is  developed.  The 
processor  is  based  on  a single-modulus  quadratic  residue 
numbering  system  where  computational  reduction  can  be 
achieved  via  an  algebraic  mapping.  Comparison  to  a 
conventional  processor  equivalent  is  made.  A 
band-selectable  DWD  for  high-frequency  spectral  analysis  is 
developed.  An  optical  guided-wave  WD  processor 
configuration  is  also  presented. 
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INTRODUCTION 


The  concept  of  a time-varying  signal,  or  one  with  a 
changing  spectrum,  has  been  given  considerable  attention. 
From  Fourier  transform  theory  it  is  known  that  frequency  is 
defined  over  an  infinite  time  interval  and  time  over  an 
infinite  frequency  interval.  Therefore,  the  concept  of  a 
joint  time-frequency  representation  of  a signal  is  not 
feasible.  Still,  a representation  of  signal  concentration 
over  time  and  frequency  serves  a useful  purpose,  an  example 
being  the  case  of  phase  modulated  signals.  As  a result, 
there  have  been  attempts  to  define  such  a representation 
based  on  the  hypothesis  of  minimum  spread.  In  1952  Page 
[1]  mathematically  defined  the  instantaneous  power  spectrum 
of  a signal.  His  representation  relied  completely  on  past 
histories  of  the  signal.  In  1963  [2]  Levin  modified  the 
work  of  Page  to  produce  an  instantaneous  spectrum  that  was 
equivalent  to  the  double  Fourier  transform  of  the  complex 
conjugate  of  the  ambiguity  function.  In  1967  Rihaczek  [3] 
combined  the  one-dimensional  energy  densities  in  the  time 
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and  frequency  domains  to  produce  a combined  time-frequency 
representation  which  he  defined  as  the  complex  energy 
density  of  the  signal.  In  1966,  Cohen[4]  developed  a 
generalized  class  of  two-dimensional  representations 
relating  electron  position  and  momentum  as  applied  to 
statistical  mechanics  then  defined  a set  of  correspondence 
rules  associating  the  marginal  distributions  of  position 
and  momentum.  The  extension  of  the  Cohen  class  to  signal 
analysis,  as  pointed  out  by  DeBruijn  [5],  is  clear.  The 
relation  between  position  and  momentum  of  an  electron,  and 
time  and  frequency  of  a signal  are  equivalent  (i.e.  they 
both  satisfy  the  Fourier  transform  relationship). 

In  this  regard,  any  member  of  the  Cohen  class  may  be 
viewed  as  an  intermediary  between  the  time  and  frequency 
domains.  Ideally,  all  properties  defining  this  class  are 
preserved  by  the  chosen  representation.  Naturally,  such  a 
case  is  not  possible  as  the  properties  are  self 

contradictory  (a  result  of  Heisenberg's  uncertainty  in 
Fourier  analysis  as  established  by  Schrodinge r ) . Still, 
most  standard  time-frequency  representations  are  members  of 
the  generalized  class  given  by  Cohen,  where  the  Cohen  class 
and  its  ideal  set  of  properties  establishes  a basis  on 
which  to  relate  time-frequency  representations  of  a signal 
[ 5-11]  . 

In  general,  most  joint  time-frequency  signal 

representations  are  obtained  via  short-time  Fourier 
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transform  (STFT)  techniques.  Fourier  transforms  are  taken 
with  respect  to  the  center  position  of  a sliding  window, 
ensemble  averaged,  and  used  to  characterize  the  spectrum. 
Such  STFT  methods  postulate  spectral  stationarity  over  the 
observation  interval.  This  defines  a compromise  between 
frequency  resolution  and  adequate  preservation  of  the 
nonstationar i ties  in  the  signal  dispersion  characteristics. 
Despite  optimum  windowing  attempts,  STFT  techniques 
generally  produce  a representation  that  is  overly  smoothed 
in  both  time  and  frequency.  Accordingly,  few  of  the  Cohen 
signal  properties  are  preserved  by  STFT  derived 
representations.  Instead,  most  properties  are  averaged 
over  the  observation  interval. 

One  signal  representation  that  does  not  necessitate 
such  a compromise  is  the  Wigner  distribution  (WD).  The  WD 
satisfies  all  of  the  properties  in  the  Cohen  class  except 
that  of  positivity.  Viewed  as  a link  between  the  time  and 
frequency  domains,  common  signal  time  domain  properties 
(instantaneous  power,  signal  energy,  instantaneous 
frequency,  and  average  frequency)  and  common  frequency 
domain  properties  (energy  spectral  density,  signal  energy, 
group  delay,  and  average  time  delay)  are  intrinsic  to  the 
WD.  The  WD  is  in  fact  the  only  member  of  the  class  of 
time-frequency  representations  to  precisely  maintain  these 
properties.  In  fact,  all  other  members  of  the  Cohen  class 
can  be  represented  by  weighted  averages  of  the  WD , 
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The  Wigner  distribution  (WD)  was  initially  developed 
by  Wigner  and  used  as  a tool  for  formulating  correction 
terms  to  the  Boltzman  transport  equations  relevant  to  low 
temperatures  [12].  Despite  the  presence  of  negative 
values,  Wigner  established  the  WD  as  an  indicator  of  the 
joint  probability  of  position  and  momentum.  DeVille  first 
applied  the  WD  to  signal  description  [13].  At  a later 
time,  DeBruijn  set  forth  a firm  mathematical  foundation  for 
the  WD[5,6].  Recently,  Claassen  and  Mecklenbrauker 
established  the  WD  and  the  discrete  WD  (DWD)  to  be  a viable 
signal  processing  tool  [7,14,15].  Their  development 
paralleled  that  of  DeBruijn,  where  a mathematical 
examination  of  the  WD  and  its  properties  formed  the  primary 
emphasis  of  their  study.  Since  that  time  several 
mathematical  descriptions  have  been  used  to  investigate  the 
WD  and  its  associated  properties  [16-23],  and  specific 
applications  of  the  WD  have  been  explored. 

Chester  examined  the  DWD  as  a speech  analysis  tool  by 
comparing  DWD  representations  of  selected  word  utterances 
with  their  associated  spectrogram  [24].  Janse  and  Kaizer 
applied  the  WD  to  loudspeakers  design  [8]  and  Preis 
examined  the  role  of  the  WD  and  associated  moments  that 
prove  important  to  audio  signal  processing  [25].  Martin 
and  Flandrin  used  a smoothed-psuedo-WD , where  a moving 
average  was  performed  on  the  WD  kernel,  to  study  biological 
rhythms  of  cave rni culous  animals  [26].  The  spWD  was  also 
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applied  to  bat  captures  by  Fandrin  and  Martin.  The 
multidimensional  WD  and  spWD  have  been  considered  in  the 
contexts  of  optics  and  image  processing  [27-29].  The  WD 
has  been  suggested  for  signal  synthesis  and  filter  design 
applications  as  well  [30,31]. 

The  increased  interest  in  the  WD  as  a signal 
processing/analysis  tool  strongly  indicates  the  need  for  a 
thorough  physical  investigation  of  the  WD  and  its  relation 
to  signal  types  in  a SP  and  DSP  setting.  Without  this 
unifying  framework  as  a basis,  applications  and  studies  of 
the  WD,  although  interesting,  are  often  inconclusive. 

One  of  the  primary  objectives  of  this  dissertation 
will  be  to  develop  such  a framework  for  the  WD.  An  in 
depth  study  of  the  WD,  and  variations  thereof,  with  respect 
to  its  role  in  the  analysis  and  description  of  signals  will 
form  the  first  phase  of  this  dissertation.  Both  DeBruijn 
[5,6]  and  later  Claassen  and  Mecklenbrauker  [7,14,15]  have 
provided  a formal  mathematical  development  of  the  WD,  and 
others  have  added  to  these  initial  studies  [7-11,16-23]. 
Chapters  2 and  3 of  this  dissertation  will  attempt  to  unify 
and  add  to  the  most  important  of  these  contributions  into 
one  comprehensive  study  as  applies  to  SP  and  DSP. 

Chapter  2 examines  where  the  Wigner  distribution 
derives  from  in  regard  to  time-frequency  representations  of 
signals  and  its  relation  to  alternative  t-f 
representations.  The  relationship  between  the  WD  as  a 
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simultaineous  t-f  signal  representation  and  the 
one-dimensional  time  and  frequency  domains  of  deterministic 
signals  forms  an  important  part  of  this  discussion.  The 
effects  of  moving  averages  performed  on  the  WD  are 
considered  in  terms  of  windowing  and  filtering  to  produce 
representations  that  maintain  weighted  averages  of  the  time 
and  frequency  domain  properties.  Window  and  filter 
selection  for  maximum  signal  concetration  are  determined. 

Chapter  3 will  apply  the  results  of  Chapter  2 to  a DSP 
setting.  The  computational  efficiency  of  FFT  data 
processing  makes  short-time  Fourier  transform  (STFT) 
techniques,  despite  their  disadvantages,  attractive  for 
real-time  implementation.  The  DWD  provides  an  alternative 
FFT  based  estimation  technique.  The  DWD  is  the  DSP 
analogue  to  the  WD  and  the  correspondence  of  the  continuous 


case  properties,  etc. 

in  a 

DSP  setting 

will 

be 

investigated.  A detailed 

comparison 

between  the 

DWD 

and 

the  conventional  FFT  based  moving  window  method,  known  as 
the  STP,  will  be  given.  The  DWD  will  be  performed  on  a 
large  class  of  signal  types  whose  dispersion  index  is  large 
and  compared  with  results  obtained  via  the  STP. 

Generally  in  DSP  applications,  a windowed  version  of 
the  WD,  termed  the  psuedo-DWD,  is  computed  resulting  in  a 
frequency  smoothed  DWD.  Flandrin[9]  offered  a more 
generalized  expression  termed  the  smoothed-psuedo-DWD 
(spDWD)  where  a weighted  average  of  the  DWD  kernel  is  used 
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to  reduce  frequency  smoothing.  The  purpose  of  the  spDWD 
was  to  effect  cross-term  reduction  in  multi-component 
signals.  In  Chapter  3 the  spDWD  is  interpreted  in  terms  of 
filtering  and  modulation  operations.  In  this  form  the 
wealth  of  knowledge  on  filtering  and  modulation  may  be 
applied  to  the  spDWD  to  improve  estimation  for  the  two 
primary  conditions  under  which  the  DWD  produces  poor 
results:  multi-component  signals  and  low  signal-to-noise 
ratios . 

Two  primary  advantages  result  from  the  spDWD:  first, 
relaxation  of  the  time/frequency  smoothing  trade-off 
associated  with  the  STP;  and,  second  the  ability  to  effect 
smoothing  similar  to  that  performed  on  the  STP  (to  improve 
performance)  without  the  computational  burden  of  computing 
multiple  FFTs  at  each  time  slice.  These  claims  are 
experimentally  tested  and  verified.  As  a final  example, 
the  spDWD  is  applied  to  estimation  of  a multipath  radio 
example  in  the  presence  of  severe  receiver  noise. 

Chapter  4 addresses  means  of  implementing  the  WD  and 
associated  variations  of  the  WD  provided  in  Chapters  2 and 
3.  Processor  configurations  tailored  to  signal  type  and 
system  application  are  provided.  General  Wigner  processing 
using  standard  FFTs  are  addressed  and  a means  of  computing 
two  N point  DWD  time  slices  via  one  N/2  point  FFT  is 
applied  to  the  case  of  real  signals.  A high-speed  DWD  (and 
spDWD)  processor  for  complex  data  is  developed.  The 
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processor  is  based  on  a single-modulus  quadratic  residue 
system  numbering  where  computational  reduction  can  be 
achieved  via  an  algebraic  mapping.  Comparison  to  a 
conventional  processor  equivalent  is  made.  A band-select 
DWD  for  high-frequency  analysis  is  developed.  An  optical 
guided-wave  WD  processor  configuration  is  presented. 

The  overall  objective  of  this  dissertation  is  to 
provide,  in  a self-contained  manner,  sufficient  preparation 
to  effectively  use  the  WD  for  signal  processing  and 
analysis  applications.  Emphasis  is  placed  on  the  DWD  (and 
spDWD)  as  a DSP  tool  and  the  theoretical  development  of  the 
DWD  is  heavily  supplemented  with  physical  examples  and 
applications  directly  related  to  signal  processing  and 
analysis.  Where  appropriate,  physical  discussions  replace 
or  supplement  formal  mathematical  descriptions. 


CHAPTER  2 

WIGNER  DISTRIBUTION — A SPECTRAL  ESTIMATOR 


2 . 1 Time-frequency  Representations 


Most  joint  time-frequency 
analysis  of  signals,  whose 
members  of  the  generalized 
representations  established  by 


representations  used  in  the 
spectra  vary  with  time,  are 
class  of  two  dimensional 
Cohen  [4]  and  defined  as 


Cf(t,(0)  = rrrf  (u+T/2)f*(u-T/2)<|>(  ^,T)e^^  ^^^dudxd^ 

JJJ  (2-1) 


where  f(t)  and  f (t)  are  the  respective  time  signal  and  its 
complex  conjugate  and  <|)(t,(A))  is  an  arbitrary  l<ernel  that 
corresponds  to  the  particular  signal  "distribution"  to  be 
considered.  The  term  "distribution"  is  set  off  in  quotes 
to  emphasize  its  loose  usage.  Assuming  C^(t,w)  is  formally 
a distribution  of  the  signal,  then  Cg(t,oo)  should  be 
positive  over  the  entire  t-w  plane  and  integration  should 
produce  the  correct  marginal  distributions.  With  respect 
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10 


to  signal  analysis  these  are 


JCj(  t,  w;  <J))dco  = I f ( t)  I ^ 
— 00 


(2-2) 


and 


00 

Jc  J ( t,co;  (f.)dt  = |F(co)  1^. 

»co 


(2-3) 


Equation  (2-2)  states  that  the  integral  of  C^(t,w)  over  all 
frequencies  at  any  fixed  time  yields  the  instantaineous 
power  at  that  time;  (2-3)  states  that  the  integral  over  all 
time  at  any  frequency  yields  the  energy  density  at  that 
frequency.  Equations  (2-2)  and  (2-3)  indicate  that 
integration  of  Cj(t,w)  over  the  whole  time-frequency  plane 
gives  the  total  energy  in  the  signal.  However,  if  C^(t,CA)) 
IS  observed  in  terms  of  signal  energy,  sharp  discrimination 
in  both  time  and  frequency  conflicts  with  the  relation 
Af~l/At  where  At  denotes  the  time  during  which  the  signal 
experiences  appreciable  change  and  Af  the  corresponding 
band  of  frequencies. 

As  an  example,  let  ( t , w) =s^ ( t , w) , where  S^(t,w) 
denotes  the  well-known  short-time  spectral  density  (STSD). 
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The  STSD  is  defined  as  the  square  modulus  of  the  short-time 
Fourier  transform  (STFT)  over  all  possible  window 
positions.  The  STFT  expresses  as 


F^(w)=Je  ^“"^f  (T)h(T-t)dT, 


(2-4) 


where  h(t)  denotes  the  moving  window  and  is  centered  about 
t.  The  power  spectrum  then  expresses  as 


S(t,co)=  |F^.(w) 


(2-5) 


Integrating  S(t,w)  over  time  at  a fixed  frequency  produces 
an  averaging  of  the  energy  density  over  the  window  as 
follows 

GO  00 

Js(t,a))dt=(l/2n)J|F(0  |^|H(w-0  l^dC.  (2-6) 

~ 00  « 00 

If,  instead,  C^(t,£A>)  is  treated  as  an  intermediary 
between  the  time  and  frequency  domains,  signal  properties 
which  characterize  the  respective  one-dimensional  time  and 
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frequency  domains  should  ideally  be  maintained  by  Cj(t,w). 
Two  of  these  properties  respectively  correspond  to  the 
instantaneous  power  and  energy  density  given  in  (2-2)  and 
(2-3).  The  remaining  properties  are  summarized  in  Table 
2-1  and  discussed  below.  Properties  2 and  6,  the  finite 
support  properties,  are 

f(t)=0  for  t<tQ  and/or  t>t^ 

^ C j ( t , co;  (() ) =0  for  t<tQ  and/or  t>tj^  (2-7) 


and 


F(w)=0  for  co<c0q  and/or  co>c0j^ 


^ C j ( t , w;  <f) ) =0  for  co<Wq  and/or 


(2-8) 


They  state  that  if  the  signal  is  time  or  bandlimited,  its 
joint  time-frequency  representation  should  also  be  time  or 
bandlimited  to  within  the  same  region.  Properties  3 and  7 
state  that  a shift  of  the  signal  in  time  or  frequency 
results  in  a corresponding  shift  in  the  time  or  frequency 
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Table  2-1. 


Characteristic 


signal  properties. 


Time 

Domain  Properties: 

1 . 

instantaneous  power 

2. 

finite  support 

3 . 

shifting 

4 . 

instantaneous  frequency  (analytic  signals) 

Frequency  Domain  Properties: 

5. 

energy  density 

6 . 

finite  support 

7 . 

shifting 

8. 

group  delay 

Signal  Distribution: 

9. 

real-valued 

10  . 

positivity 

14 


of  Cj(t,w)  as  follows; 


g(t)  = f(t-tQ)  =>  Cg(t,o);<f.)=C^(t-tQ,w;<|))  ; (2-9) 


and 


G(co)=F(  w-cOq)  Cg(t,oo;4.)=Cf  (t,w-a)Q;(())  . (2-10) 


Property  4 states  that  the  average  frequency  of  C^(t,w) 
any  time  is  equal  to  the  derivative  of  the  phase  of 
signal,  or 


J‘ooC£(  t,  (o;  <|))dco 

!! = -d0(t)/dt,  (2-11) 

00 

JCj(  t,w;<f))da) 


where  the  envelope  of  f(t)  is  of  the  form  a(t)e^®^^^. 
real  signals,  the  average  frequency  is  zero  and 
analytic  signals,  the  average  frequency  is  equal  to 
instantaneous  frequency.  Property  8 states  that 
average  time  of  Cj(t,w)  is  equal  to  the  group  delay. 


at 

the 


For 

for 

the 

the 


That 
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is,  if  Tg{  co)=-d(|)(  co)/dw,  where  F(  co)  =R(  co) e^  ^ , then 


J‘tCj(t,w;(j))dt 

” = T (CO).  (2-12) 

OD  y 

JCj(  t,co;4>)dt 


Property  9, 

Cj(  t,w;(j))  = C*  j(  t,co;  (f>)  , 


(2-13) 


states  that  C^(t,co)  is  real-valued.  Property  10, 

Cf  (t,o);<f))>0,  (2-14) 

states  that  C^(t,oo)  is  everywhere  positive. 

Although  ideal,  a generalized  time-frequency 

representation  for  which  all  the  properties  in  Table  2-1 
hold,  reguardless  of  the  signal,  is  not  possible. 

Consider,  for  example,  a representation  that  satisfies 
properties  2 and  10:  then  from  (2-7),  if  the  signal  is  0 

for  t<0,  Cj(t,co)=0  for  t<0.  By  definition,  the  time 
average  of  C^(t,co),  where  Property  10  defines  C^(t,co)  as 
everwhere  >0,  is  >0  at  all  frequencies.  Obviously, 
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Property  8 (Equation  2-12)  in  this  case  does  not  hold  as 
the  group  delay  for  causal  signals  is  not  strictly 
positive.  Although  these  properties  may  not  absolutely 
hold,  they  may  be  satisfactorily  maintained  under  certain 
signal  constraints. 

Again,  consider  the  case  where  C£(t,co)  corresponds  to 
the  STSD.  Integrating  the  STSD  over  time  at  an  arbitrarily 

fixed  frequency  was  shown  to  produce  a windowed  average  of 
the  energy  density  of  the  form  (2-6).  In  suit,  the  first 
moment  over  time  of  S(t,(o)  (where  for  simplicity  the 
subscript  has  been  omitted)  produces  an  averaging  of  the 
group  delay.  That  is 

00  CO 

;ts(t,co)dt  ;t  (OF'(w,odc 

— GO  — 00  ^ 

CO  00 

JS(t,w)dt  ;FMw,OdC 

— 00  —00 

f 

where 

FMoo,a  = |f(j:)  i^|H(co-a  i^.  (2-15) 

If  the  significant  spectral  contributions  of  the  signal  do 
not  change  appreciably  over  the  length  of  the  window,  then 
(2-15)  and  (2-16)  correspondingly  do  not  vary  greatly  over 
the  window  length  and  the  respective  average  values  are 
sufficiently  close  to  the  instantaneous  values.  A signal 


17 


of  this  type  is  often  referred  to  as  quasi-stationary 
or  stationary  over  the  window.  Conversely,  a signal  whose 
dispersion  characteristics  undergo  significant  change  over 
the  window  length  is  referred  to  as  a time-varying  or 
r^onstationary  signal.  Naturally,  the  STSD  would  not 
be  an  adequate  representation  for  signals  that  are 
time-varying  and  an  alternative  choice  for  C^(t,w)  must  be 
considered.  More  specifically,  the  objective  in  selecting 
C£(t,w)  is  to  optimally  represent  signal  concentration  over 
time  and  frequency,  which  is,  in  fact,  limited  by  the 
uncertainty  relation  of  Heisenberge,  and  the  Cohen  class 
with  its  ideal  set  of  properties  serve  as  a viable  basis  on 
which  to  relate  time-frequency  representations  of  a signal 
[6,7,8,11]  . 


2 . 2 Wigner  Distribution 

The  continuous  Wigner  distribution  (WD)  of  a signal 
f ( t ) is  defined  as 


00 

W(t,w)  = Jf  (t+T/2)f*(t-T/2)e~^“'^dT  (2-16) 

_ 00 


and  is  obtained  from  (2-1)  by  setting  <j>=l . The  WD 
satisfies  all  of  the  defining  properties  in  Table  2-1 
except  that  of  positivity  [6].  (Note  should  be  made  that 
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Property  6 is  constrained  to  a single  bandpass  region;  the 
case  of  multiple  nonzero  bands  is  discussed  in  Section 
3.3.)  The  occurrence  of  negative  values  that  allows  the 
other  nine  properties  to  simultaneously  hold,  however, 
prevents  the  WD  fron  being  described  formally  as  a 
distribution.  The  existence  of  the  negative  values  may  be 
equated  to  that  of  negative  energy  states  introduced  in 
quantum  mechanics  or  of  negative  temperatures  in 
thermodynamics . 

Hudson  [23]  has  proven  the  WD  as  applied  to  quantum 
mechanics,  in  the  form  of  a phase-space  distribution,  to  be 
everywhere  positive  if  and  only  if  the  Schrodinger  state 
vector  is  the  exponent  of  a quadratic.  With  respect  to 
signal  estimation,  Hudson's  proof  equates  to  positivity  for 
the  WD  given  the  signal  envelope  is  the  exponent  of  a 
quadratic.  In  an  earlier  paper,  DeBruijn  proved 
application  of  constrained  Gaussian  smoothing  operators  on 
the  WD  could  effectively  elliminate  the  negative  value 
problem  [ 5 ] . 

In  addition  to  those  properties  which  characterize  the 
Cohen  class,  the  WD  has  a number  of  desirable  attributes 
that  make  it  well  suited  to  signal/image  processing 
applications  [6,14,15].  For  example,  two  common  linear 
operations  on  signals  are  filtering  and  modulation.  The 
duality  of  each  operation  with  respect  to  time  and 
frequency  is  evident  in  the 
specifically. 


resulting 


WD. 


More 
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filtering:  (2-17) 

00 

y(t)  = x(t)*^h(t)  Wy(t,w)  = Jw^(t-T,w)Wj^(t,oo)dT 

_ 00 

modulation:  (2-18) 

00 

y(t)  = x(t)h(t)  =>  Wy(t,o))  = Jw^(  t,w)Wj^(  t,  w-v)dv 

— 00 

where  W^,  , Vl^,  are  the  respective  WD's  of  y,  x,  h and 

denotes  convolution.  These  two  operations  are  important  to 
real-time  processing  of  the  WD.  Another  characteristic  of 
the  Wigner  is 


Wj(t,w)  = Wp,(w,-t)  , 


(2-19) 


where  is  the  WD  of  the  signal  f(t),  as  defined  in 
Equation  (2-16)  and  is  the  WD  of  the  signal  transform 
with  integration  over  frequency  as  follows, 


00 

Wp(w,t)  = jF(w+^/2)F*(w-£;/2)e“^"‘^d^,  (2-20) 

— CO 


where  F(co)  denotes  the  Fourier  transform  of  f(t). 


The 


20 


duality  between  t and  w is  evident  in  (2-19).  This 
demonstrates  again  the  validity  of  the  WD  as  a 
time-frequency  representation. 

2.2.1  Windowed  Wigner  Distribution 

In  practical  applications,  a windowed  version  of  the 
WD  is  used.  If  the  window  is  of  an  appropriate  Gaussian 
form  the  WD  becomes  observably  nonnegative  (Section  2.2.3). 
The  WD  of  a windowed  version  of  the  signal  follows  from 
(2-18).  If  h(t)  denotes  a window,  centered  at  t',  applied 
to  the  signal,  f(t),  such  that  y ( t ) =f ( t ) h ( t-t ' ) , then 


The  WD  can  be  seen  as  a function  of  the  window  position. 
Sliding  the  window  along  the  time  axis  yields  a series  of 
WDs.  The  pseudo-WD , or  pWD,  is  defined  by  observing 
the  WD  only  at  the  window  center  t=t'.  As  the  window  moves 
along  the  time  axis,  it  is  expressed  as 


CO 


(2-21) 


GO 


(2-22) 


— CD 
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or,  equivalently, 


W 


(t,w)  = [f  (t+x/2)f*(  t-T/2)h(T/2)h*(-T/2)e 

(2-23) 


The  pWD  produces  spreading,  or  smoothing,  in  frequency 
that  is  equivalent  to  a convolution  in  the  frequency  domain 
of  the  WD  of  the  signal  and  the  WD  (at  t=0)  of  the  window. 
Although  the  time  domain  properties  in  Table  2-1  are 
unaffected,  the  frequency  smoothing  translates  to  an 
averaging  over  the  WD  of  the  window  for  properties  5,  6 and 
8.  For  practical  applications  there  is  no  loss  of 
generality  in  assuming  the  window  function  to  be 
real-valued,  in  which  case  Property  4 will  be  maintained 
precisely  by  the  pWD . For  a causal  system  with  impulse 
response  approaching  zero  outside  a finite  time,  Tq,  (as 
for  example  loudspeakers  [8])  the  pWD  closely  resembles  the 
WD  if  the  window  length  exceeds  Tg  (i.e,  properties  5,6 
and  8 almost  hold) . 

2.2.2  Doubly  Weighted  Wigner  Distribution 


Frequency  spreading  of  the  pWD  relative  to  each  time 
slice  can  be  reduced  via  a moving  filter.  That  is,  using 
the  t-w  duality  of  the  WD  as  given  in  (2-19),  each  filtered 


pWD  slice  is  satisfied  by  (2-17)  where  the  convolution 
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function  Wj^  becomes  the  WD  of  the  filter  at  time  t=0.  In 
this  form  the  spWD  can  be  treated  as  a doubly  biased(or 
smoothed)  WD  where  the  relative  frequency  and  time 
resolutions  can  be  controlled  by  the  choice  of  window  and 
filter.  This  is  the  one-dimensional  counterpart  to  the 
composite  WD  in  [27]. 

To  compute  the  weighted  WD,  the  signal  is  windowed  over 
a finite  time  centered  about  an  arbitrarily  fixed  time  tQ, 
producing  f ( t ) = f ( t ) h^^  ( t-tg  ) . The  transform  of  the  windowed 
function  is  computed  then  multiplied  by  a filter  centered 
about  a frequency  Wq , to  yield 

^t0,(o0^“^  = { F(  03)  ( w)  exp(  - jcotQ  ) }H2  ( oo-Wq  ) (2-24) 


or,  equivalently. 


^t0,o)0^^^  = { f ( t )h^  ( t-tg  ) } *h2  ( t ) exp(  josgt ) . (2-25) 


The  weighted  WD  is  defined  by  observing  the  WD  at  t = 
03=0)q  and  at  each  point  (t,o>)  expresses  as 


W 


^(t,o))  = ( 1/2 II ) JJWj  ( z , V ) Wj^^  ( 0 , 03-v  ) Wj^2  ( t-z  , 0 ) dzdv  . 

(2-26) 
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As  the  equivalent  of  windowing  and  filtering  is 
performed  on  the  signal  to  produce  the  wWD,  a weighted 
average  of  both  the  time  and  frequency  domain  signal 
properties  is  maintained  by  the  wWD.  In  fact,  all  other 
members  of  the  Cohen  class  can  similarily  be  represented  by 
weighted  averages  of  the  WD  in  the  following  manner  [6,14], 


0000 

Cf(t,w;<^)  = (l/2ii)JJWj(C,^)e(t-C,w-^)d<;d^,  (2-27) 

CD  00 


for  W^(t,w)  the  WD  of  f(t)  and  0 the  appropriate  weighting 
where 


0000 

9(t,co)  = 

0000 


(2-28) 


for  (J)(t,w)  satisfying  (2-1).  For  example,  setting 
9(  t , co)  =e^  ^ in  (2-27)  would  yield  the  double  FT  of 
W^(t,w),  which  would  in  fact  produce  the  ambiguity  function 
of  the  signal  (Section  2.3). 
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2.2.3  Maximum  Concentration  Weighting  Functions 

A natural  question  arises  as  to  the  optimum  selection 
of  window  and  filter  functions.  As  known,  signal 
concentration  in  the  t-w  plane  is  limited  by  Heisenberg's 
uncertainty  inequality  [5] 


CO  00 

[ Jt^|f(t)  |^dt]^/^[  Ja)^F(w)dco]^/^>|  |f  I |/4n 

— 00  _oo 


(2-29) 


where  | | ' | | denotes  the  L2  norm  and  the  lobe  center  is 
assumed  to  be  at  the  origin.  For  the  case  of  the  WD 
(properties  1 and  5)  Equation  (2-29)  translates  to 

00  GO 

[t^W(  t,  co)dtdcA)rco^W(  t,co)dtdw>(  | |f  | |/4n)^ 

(2-30) 


Both  second  global  moments  are  strictly  positive  for  the  WD 
[14]  and,  therefore,  (2-30)  can  be  interpreted  as  the 
product  of  the  variances  of  W(t,w)  about  the  lobe  center 
over  time  and  over  frequency.  A well  known  fact  from 
quantum  mechanics  is  that  the  inequality  in  (2-29)  becomes 
an  equality  if  the  signal  waveform  is  a Gaussian  (i.e.  the 
measure  of  uncertainty  is  a minimum).  Furthermore, 
Heisenberg's  inequality  approaches  an  equality  if  the 
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waveform  uniformly  approaches  a function  for  which  it 
becomes  an  equality  [5].  Therefore,  a logical  conclusion 
would  be  to  select  the  window  to  approach  a Gaussian  shape. 
From  (2-19)  and  (2-20)  the  same  conclusion  holds  for  the 
selection  of  the  moving  filter. 

From  (2-26)  the  weighting  function  (or  moving  average) 
in  (2-27)  for  the  weighted  WD  is  a separable  function  of 
the  form  0^(t)02(w).  From  (2-16)  and  (2—26),  if  the  window 
is  a Gaussian  with  variance  3^,  the  corresponding  smoothing 
function,  0-(w),  is  also  a Gaussian,  with  variance  2/6  ; 
likewise,  from  (2-19),  (2-20)  and  (2-26),  if  the  filter  is 
a Gaussian  with  variance  6j,  the  smoothing  function  0^^  ( t ) 
is  a Gaussian  with  variance  2/6£.  For  the  case  of  Gaussian 
windowing  and  Gaussian  filtering  on  the  signal  then  the  wWD 
satisfies  (2-27)  where 


0(  t,  w)  = it(  6^3f 


^/^exp[-6ft^/4-6^w^/4] 


(2-31) 


which  is  the  Weierstrass  transform  and  converges  rapidly. 
Furthermore,  if  the  window  and  filter  are  selected  such 
that  6^3j<16,  then  the  corresponding  smoothing  functions 
satisfy  the  positivity  requirement  established  by  DeBruijn 
[5]  and  the  weighted  WD  becomes  everywhere  positive. 
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2 . 3 Ambiguity  Function 

The  double  Fourier  transform  of  the  WD  is  the 
ambiguity  function  (AF),  which  is  common  to  radar  and  sonar 
applications.  The  AF  expresses  as 

00 

A(t,w)  = Jf (T+t/2)f*(T-t/2)e“^“^dT  (2-32) 

— 00 

and  can  be  obtained  from  (2-1)  by  setting  <(>=2  n5  ( r-t ) 6 ( Jl-w)  . 

Certain  properties  of  the  AF  prohibit  ideal  resolution 

between  targets.  Specifically,  these  properties  require 

that  the  maximum  height  and  the  total  volume  of  the  AF  are 

2 

both  equal  to  (2E)  where  E is  the  energy  in  the 
doppler-shif ted  signal.  In  effect  a finite  height  with 
finite  volume  is  used  to  represent  the  ambiguity  diagram 
possibly  preventing  resolution  between  two  closely  spaced 
targets.  The  width  of  the  peak  may  be  narrowed;  however, 
the  height  of  the  AF  must  be  raised  elsewhere  to  retain  the 
(2E)  volume  causing  unwanted  peaks  (pillow  effect)  [33]. 
As  such  a trade-off  between  doppler  frequency  resolution 
and  range  ambiguity  must  be  made. 

The  ambiguity  function  has  been  applied  extensively  to 
both  radar  and  sonar.  Its  plot,  called  the  ambiguity 
diagram,  serves  to  represent  the  response  of  the  matched 
filter  to  the  signal  to  which  it  is  matched  and  to 
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doppler-f requency-shif ted  signals.  As  such,  the  AF  can  be 
used  as  an  indicator  of  the  limitations  and  advantages  of 
particular  waveform  classes.  The  AF  is  a representation  of 
time  delay  and  frequency  shift  and  its  properties  should  be 
interpreted  accordingly.  Consider,  for  example,  the 
shifting  properties  in  Table  2-1.  Time  delay  and  frequency 
translation  of  the  signal  directly  correspond  to  respective 
time  and  frequency  shifts  in  the  WD.  In  the  AF,  time  delay 
of  the  signal  manifests  as  a translation  in  frequency  and 
frequency  translation  appears  in  the  AF  plane  as  a shift  in 
time.  This  clearly  agrees  with  the  double  Fourier 
transform  relationship  between  the  WD  and  the  AF . 

The  Fourier  transform  relationship  between  the  AF  and 
the  WD  can  be  used  to  translate  the  wealth  of  knowledge 
associated  with  the  AF  to  the  WD,  keeping  in  mind  the 
transform  relationship.  Consider  the  following  example,  a 
Gaussian  waveform  is  known  to  produce  maximum  ambiguity  (or 
variance)  in  the  AF  plane;  the  double  Fourier  transform 
relationship,  then,  would  indicate  the  Gaussian  to  produce 
minimum  variance  for  the  WD,  which,  in  fact,  is  the  case. 

2 . 4 Short-time  Spectral  Density 

The  short-time  spectral  density  (STSD),  defined  in 
(2-4)  and  (2-5),  can  be  obtained  from  (2-1)  by  setting  <j> 
equal  to  the  AF  of  the  window,  h(t),  and  satisfies  only 
properties  3,7,9  and  10.  As  discussed  in  Section  2.1,  if 
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the  signal  properties  do  not  vary  greatly  over  the  interval 
in  which  the  spectrogram  is  being  taken,  the  remaining 
properties  are  suitably  averaged ; this  is  generally  true 
for  short-time  stationary  signals  where  the  signal  is 
assumed  stationary  over  short  time  intervals  and  the  STFT 
is  computed  over  the  interval  then  concatenated  in  time  to 
produce  the  spectrum.  For  signals  characterized  by 
instationarities  in  the  signal,  which  is  true  for  audio 
signals  (hence  speech)  [24],  this  quasi-stationary 
assumtion  can  mask  important  information  contained  in  the 
signal . 

A comparison  can  now  be  made  between  the  spectrogram 
expressed  in  this  manner  and  the  smoothed-psuedo-WD ; for 
the  spectrogram  smoothing  in  both  the  time  and  frequency 
domains  are  interdependent  leading  to  the  well-known 
problem  of  trade-off  between  spectral  resolution  and  the 
quasi-stationary  assumption  of  the  signal.  In  the 
weightedWD  this  problem  is  avoided  by  allowing  time  and 
frequency  smoothing  to  be  computed  independently,  that  is 
9(t,w)  in  (2-27)  is  a separable  function  of  the  form 
0^(t)02(w)  for  the  wWD. 


Another 

common 

two-dimensional 

representation 

for 

audio  signals  is 

the  cumulative 

spectra 

(CS)  . 

The 

cummulative 

spectra 

is  actually  a 

special 

type 

of 

spectrogram 

where 

the  only  proper 

ties  in 

Table 

2-1 

satisfied  by 

the  CS 

are  also  3,7,9  and 

10. 

CHAPTER  3 

APPLICATIONS  TO  DIGITAL  SIGNAL  PROCESSING 


3 . 1 Discrete  Wigner  Distribution 

Generally  estimation  of  a time-varying  spectra  is 
based  on  the  short-time  Fourier  transform  (STFT).  Here  the 
process  is  assumed  stationary  over  short  time  intervals  and 
the  discrete  Fourier  transforms  (DFTs)  of  these  intervals 
are  ensembled  then  used  to  characterize  the  spectrum.  STFT 
methods  postulate  stationarity  over  the  observation 
interval;  this  defines  a compromise  between  frequency 
resolution  and  adequate  preservation  of  the 
nonstationarities  in  the  signal  dispersion  characteristics. 
Despite  optimum  windowing  attempts,  STFT  techniques 
generally  produce  a representation  that  is  overly  smoothed 
in  both  time  and  frequency.  Still  the  computational 
efficiency  of  FFT  data  processing  makes  STFT  techniques 
attractive  for  real-time  implementation  and  the  short-time 
periodogram  (STP),  defined  as  the  square  modulus  of  the 
STFT,  is  a common  tool  used  to  characterize  time-varying 
spectra.  The  discrete  form  of  the  WD,  the  DWD,  uses  an  FFT 
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and  serves  as  a 
methods  (namely  the 
In  this  section  the 
is  a competative 
comparison  between 
fully  developed. 

The  discrete 
defined  as 


viable  alternative  to  conventional  STFT 
STP  and  modified  (or  smoothed ) -STP ) . 
DWD  will  be  developed.  Because  the  DWD 
alternative  to  the  STP,  a detailed 
each  will  be  made  after  the  DWD  has  been 

time  Wigner  distribution  ( DTWD ) is 


W(n,w)=2r  X ( n+m) X* ( n-m) e 
m=-oo 


(3-1) 


where,  as  is  common  practice,  n designates  nT,  for  T the 
time  between  samples.  By  observing  the  "power  of  two"  of 
the  exponent,  the  DTWD  can  be  seen  to  be  periodic  in 
frequency  with  period  n/T  (whereas  the  Fourier  spectrum  is 
periodic  with  period  2n/T)  which  can  lead  to  aliasing  of 
the  DTWD  in  the  frequency  domain.  An  alternative 
definition  is  possible  which  omits  the  power  of  two  in  the 
exponent.  However,  by  doing  this  one  must  recall  that 
frequencies  at  w in  the  signal  are  translated  out  to  nearly 
2(A)  during  computation  of  the  inner  product  function,  or  WD 
kernel,  x ( n+m ) x* ( n-m ) . For  clarity,  consider  a sinusoid 
with  frequency  fQ=375Hz  (Figure  3-1):  in  the  inner  product 
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Figure  3-1.  STFT  of  a single  sinusoid  at 
fQ=375Hz  and  sampled  at  IkHZ. 


Figure  3-2. 


DWD  of  a single  sinusoid  at 

f„=375Hz  sampled  at 

iRHz  where  the  tone  at  2wq 

is  an  aliased  spectral  component. 
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Figure  3-3.  DWD  of  the  375Hz  sinusoid  sampled 
at  2kHz. 
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£q  is  mapped  to  dc  and  2f^,  which  folds  into  the  t-f  plane 
interval  as  an  aliased  contribution  (Figure  3-2).  (Notice 
dc  information  is  difficult  to  obtain  from  the  WD  due  to 
the  spectral  "corruption"  around  dc  associated  with  the 
inner  product. ) To  prevent  aliasing  in  the  frequency 
domain,  the  signal  must  be  sampled  at  a rate  commensurate 
with  the  effective  bandwidth  of  the  WD  kernel.  This  can  be 
guaranteed  by  sampling  the  signal  at  double  the  Nyquist 
rate . 


By  defining  the  DTWD  as  in  (3-1),  frequencies  that 
occur  at  w in  the  signal  also  occur  at  w in  the  DTWD.  That 
is,  frequencies  translated  out  to  2co  during  computation  of 
the  inner  product  are  isomorphically  folded  back  to  w 
(Figure  3-3).  If  the  signal  has  been  sampled  at  twice  the 
Nyquist  rate  the  DTWD  for  each  period  is  zero  over  the 
second  half  of  the  interval,  n/T<co<2 n/T ; hence  no  frequency 
aliasing  occurs.  If  supplied  with  a discrete  time  sequence 
sampled  at  the  Nyquist  rate,  interpolation  of  the  data  can 
be  used  to  satisfy  this  criterion.  Observe  from  Figures 
3-2  and  3-3,  the  DWD  introduces  spectral  components  around 
dc.  Where  this  inherent  dc  "corruption"  derives  from  and 
its  consequential  elimination  are  discussed  in  Section  3.3. 

Alternative  definitions  for  the  DTWD  have  been 
considered  [18,22,34],  each  of  particular  merit;  however, 
the  need  to  sample  at  twice  the  Nyquist  rate  to  combat  the 
problem  of  aliasing  has  been  common  to  all.  That  is,  the 
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DTWD  of  the  oversampled  signal  is  fully  described  over 
(0,n/T).  Using  the  definition  in  (3-1)  and  assuming  the 
signals  to  be  sampled  at  twice  the  Nyquist  rate,  the  Cohen 
class  properties  in  Table  2-1  can  be  applied  to  the  DTWD. 
Properties  1 and  5,  the  instantaneous  power  and  energy 
density,  defined  by  (2-2)  and  (2-3),  respectively,  in 
Chapter  2,  express  as 

00 

Jw(n,w)dw  = |x(n) (3-2) 

~00 


and 


Z W(n,w)  = ( |F(w)  |^+|F(«+il/T)  I )^,  (3-3) 

n=-oo 

where  the  second  term  in  (3-3)  vanishes  for  the  cases  of 
analytic  or  oversampled  signals.  The  finite  support 

property  with  respect  to  time  (2)  follows  in  a 
straightforward  manner  from  (2-7)  and  expresses  as 

x(n)=0  n>nQ  and  n<n^^ 

[W(n,co)  ]modN=0,  n>np  and  n<n^ . 


( 3-4) 
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The  finite  support  property  with  respect  to  frequency 
(Property  6,  Equation  (2-8)  is  not  as  straightforward:  for 

analytic  signals,  where  the  spectrum  vanishes  between  n/T 
and  2 n/T , the  following  relation  holds, 

X(w)  = 0 co^<co<n/T 

=>[W(  n,  co)  ]modn/T=0  C0j^<n/T;  (3-5) 


and,  for  oversampled  signals,  where  the  spectrum  vanishes 
between  n/2T  and  3n/2T,  the  following  relation  is 
satisfied, 

X((o)=0  a)^<(A)<W2,  “i<^/2T,  W2>3n/2T 

=>W(n, (A))modn/T=0  (o^<(o<oo2-n/T.  (3-6) 


As  in  the  continuous  Wigner,  (3-5)  and  (3-6)  exclude  the 
case  of  multiple  bands,  which  will  be  considered  separately 
in  Section  3.3.  The  shifting  properties  3 and  7 (Equations 

(2-9)  and  (2-10))  correspond  directly  to  the  DTWD 
respectively  as  follows: 


x(n-nQ)=^W(n-nQ,w)  ; 


(3-7) 
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and 


X(w-(»)Q)=>W(n,w-coQ)  . 


(3-8) 


The  group  delay  of  the  di screte-time  sequence  (Property  4) 
follows  analogously  from  (2-12)  as  below 


The  notion  of  an  instantaneous  frequency  (Property  8) 
obviously  cannot  apply  to  discrete-time  sequences;  however, 
the  instantaneous  frequency  defined  for  cont inuous-time 
sequences  and  satisfying  (2-11)  can  be  translated  to  the 
case  of  the  DTWD  by  observing  the  forward  and  backward 
differences  of  the  phase  about  n [15].  Here,  the  following 
equality  holds 


N-1 

Z nW  ( n , 00 ) 


N-1 

Z W(n,oo) 
n=0 


(3-9) 


] modit/T 


n/2T 

J ooW(n,w)doo 


I 


JI/2T 

J W(  n , oo)  doo 


-II/2T 


( 3-10) 
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The  real-valued  property  (9)  of  Equation  (2-13)  expresses 
for  the  DTWD  as  below. 


W(n,w)=^W  (n,w) 


(3-11) 


As  with  the  continuous  case,  property  10  (Equation  (2-14)), 
that  of  positivity  does  not  hold  for  the  DTWD. 

For  finite  duration  sequences,  of  length  NT,  the  DTWD 
can  uniquely  be  expressed  by  its  samples  taken  at  intervals 
of  (0  = kix/NT  as  below. 


N/2 

W(n,k)  = 2Z  x(n+m)x*(n-m)e'^"'^^'^/^, 
m=-N/2+l 


3-12) 


which  is  periodic  in  w=n/T.  For  analytic  signals  where  the 
even  or  odd  samples  fully  describe  the  signal,  hence  making 
the  spectrum  nonzero  over  only  one-half  the  interval  (i.e. 
single-sideband  SSB),  the  signal  need  only  be  sampled  at 
the  Nyquist  rate  to  prevent  aliasing. 

At  this  stage  an  important  point  needs  to  be  made. 
Although  aliasing  in  the  frequency  domain  can  be  prevented 
by  sampling  the  signal  at  a rate  commensurate  with  the 
effective  bandwidth  of  the  Wigner  kernel,  aliasing  in  the 
DWD  plane  is  unavoidable.  This  fact  becomes  clear  by 
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recalling  the  inverse  relation  Af~l/Zit  discussed  at  the 
beginning  of  Chapter  2.  A simple  procedure  to  minimize 
aliasing  while  maintaining  the  relative  simplicity  of 
computing  the  DWD  is  to  oversample  the  signal  as  stated, 
window  the  signal  with  a Gaussian  shaped  window  to  minimize 
the  signal  BT  product  (this  will  make  the  DWD  observably 
nonnegative  as  well)  and  zero  pad  the  window  tails  to 
reduce  aliasing  in  time  (or,  equivalently,  select  the 
window  variance  to  be  relatively  small).  For  a fixed 
transform  length,  windowing  and  zero  padding  will  reduce 
frequency  resolution  at  each  time  slice.  Recall,  however, 
for  the  DWD  increasing  the  window  length  to  improve 
frequency  resolution  does  not  degrade  the  temporal 
characteristics  of  the  DWD.  Therefore,  the  amount  of  zero 
padding  is  limited  primarily  by  the  processor  capabilities. 


3 . 2 Windowed  Discrete  Wigner  Distribution 

For  real-time  applications  data  sequences  are 
frequently  long  in  which  case  the  pseudo-DWD  is  computed  as 

W(n,k)=2i:  x(n+m)x*(n-m)  |h(m)  I (3-13) 

m=-L+l 

where  h(m)  is  a window  of  length  2L-1  centered  about  n. 
The  pseudo-DWD  is  usually  called  simply  the  DWD  and,  unless 
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stated  otherwise,  DWD  will  refer  to  the  pseudo-DWD  of 
(3-13).  The  pDWD  is  defined  analogous  to  its  continuous 
domain  counterpart,  the  pWD  and  is,  by  definition,  the 
concatenation  of  DWD  slices  over  time  where  each  slice, 
taken  with  respect  to  the  window  center,  expresses  as 

W (0,k)=22:^x(m)x*(-m)  |h(m)  | ^exp“^  ^ (3-14) 
m=-L+l 


and  n designates  the  absolute  center  position  of  the  moving 
window.  Each  DWD  slice  given  above  can  be  shown  to  be  the 
DFT  of  a rearranged  inner  product  function,  ij^'^(m), 
frequently  referred  to  as  the  DWD  kernel.  The  DWD  kernel 
is  obtained  by  rearranging  the  inner  product, 
x(m)h(m)x  (-m)h  (-m),  such  that  terms  in  (3-14)  from  zero 
to  N/2  form  the  first  half  of  the  DWD  kernel  and  the  terms 
in  (3-14)  from  -N/2+1  to  -1  form  the  second  half  of  the 
kernel  and  i^  (m)  is  defined  to  run  from  0 to  N-1.  Using 
the  relation  exp ( - jmk2 n/N ) = exp ( - j ( m+N ) k2 n/N ) , the  DWD  can 
now  be  expressed  as 

W(n,k)=2E  i^’^(m)exp“^™'^^'^/^ 
m=0 


(3-15) 
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where  i^’^(m)  designates  the  DWD  kernel  as  defined  above 
taken  with  respect  to  a window  centered  about  the  point  n. 
Expressed  in  the  form  of  (3-15),  the  DWD  is  just  the  DFT  of 
the  kernel  i^'’(m),  hence  lending  itself  readily  to  current 
FFT  implementations. 

Due  to  the  symmetry  of  the  DWD  kernel,  data 
rearrangement  is  not  actually  needed  in  the  DWD  computation 
if  the  signal  is  real.  That  is,  the  DFT  of  the  windowed 
sequence  by  definition  assumes  periodicity  of  the  signal 
outside  the  observation  interval  and  application  of  the  DWD 
kernel  directly  to  the  DFT,  minus  data  rearrangement,  is 
equivalent  to  starting  the  period  at  a different  point  in 
the  periodic  sequence.  For  practical  applications  there  is 
no  loss  of  generality  in  assuming  the  window  to  be  a real 
symmetric  function  in  which  case  the  DWD  expresses  as 


N-1 

W(n,k)=2Z  |h(m) 
m=0 


,„^„„„-jmk2n/N 
( m ) exp 


(3-16) 


where 

rectangular 


is  the  unrearranged  DWD  kernel  given  a 
window  of  length  N and  is  by  definition 


symmetric  about  n. 
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3.2.1  Comparison  with  the  Short-time  Periodogram 

The  DWD  may  now  be  compared  directly  with  the  popular 
short-time  periodogram  (STP).  The  STP,  defined  as  the 
square  modulus  of  the  short-time  Fourier  transform, 
expresses  as 


P(n,k)=l/Nj  L x(m)h(n-m)e  j2fimk/N|2^  (3-17) 

' m=0  I 

where  h(n)  is  a window  of  length  N.  Despite  its 
limitations,  the  STP  can  easily  be  computed  using  FFTs  and, 
as  such,  has  been  a standard  estimation  procedure.  The 
primary  limitation  of  the  STP  has  been  the  trade-off 
between  frequency  resolution  and  preservation  of  signal 
instationarities  in  selecting  the  appropriate  transform 
length.  Consider  for  example  a square  pulse  envelope  on  a 
400Hz  carrier  as  illustrated  in  Figure  3-4:  in  the  STP 

(Figures  3-5  and  3-6)  the  lobe  width  increases  at  the  edges 
where  the  quasi-stationary  assumption  over  the  window 
becomes  invalid.  For  the  longer  transform  length  frequency 
smoothing  is  reduced,  as  evident  by  the  narrowing  of  the 
spectral  lobes;  however,  the  transient  part  of  the  signal 
becomes  significantly  degraded.  Application  of  a window  as 
shown  for  the  case  of  the  Hamming  window  (Figures  3-7  and 
3-8)  can  be  used  to  improve  the  spectral  output;  however. 


VALUES 


2tM8 


Figure  3-4. 


Square  pulse  envelope  on  a 
400Hz  carrier. 
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Figure  3-5.  STP  of  square  pulse  modulation 

example  for  a 64  point  transform 
length . 


Figure  3-6 . 


STP  of  square  pulse  modulation 
example  for  a 256  point  transform 
length . 
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Figure  3-7.  STP  of  square  pulse  modulation 

example  for  a 64  point  transform 
length  with  Hamming  windowing. 
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STP  of  square  pulse  modulation 
example  for  a 256  point  transform 
length  with  Hamming  windowing. 


Figure  3-8. 


45 


improvements  are 
changes . 

The  inherent 
by  observing  the 
discussed  in  Chapt 
<fr(t,w)  in  (2-1) 
STSD  which  assumed 
the  double  Fouri 
( Equation  ( 2-28 ) ) 
can  be  expressed 
DWD  of  (3-12)  whe 
nonseparable  funct 


limited  primarily  to  minor  temporal 


t-co  trade-off  limitation  becomes  apparent 
STP  as  the  DSP  analog  to  the  STSD 
er  2.  Recall,  the  weighting  function 
was  a nonseparable  function  for  the 
the  form  of  the  AF  of  the  window.  Using 
er  transform  relation  between  ()>  and  9 
and  between  the  AF  and  the  WD,  the  STP 
as  a weighted  average  of  the  (unwindowed) 
re  the  weighting  term,  0(n,k),  is  a 
ion  of  the  form 


N-1 

0„(n,l<)  = Z h(n+m)h(n-m) 
^ m=0 


g- j2jimk/N 


( 3-18) 


For  the  pDWD  the 
windowing  assumes  the 


weighting  function 
following  form 


resulting  from 


N-1 

9^(k)=E  h(m)h(-m 
m=0 


- j2nmk/N 

® 9 


(3-19) 


which  is  the  DWD  of  the  window  at 
(2-22)).  Here  the  weighting  term  is 


time  zero  (Equation 
a function  of  k only. 
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where  windowing  of  the  DWD  does  not  smooth  the  time  domain 
signal  properties  avoiding  the  standard  time/frequency 
trade-off  problem  associated  with  the  STP.  Figures  3-9  and 
3-10  give  the  DWD  for  respective  64  and  256  point 
transforms  of  the  square  pulse  modulation  in  Figure  3-4. 
In  the  DWD  the  lobe  widths  remain  relatively  constant  over 
the  time  interval.  As  with  the  STP,  appropriate  windowing 


of  the  data  can  be 

used 

to 

improve 

the  estimation 

particularly  for  the 

case 

of 

Gaussian 

smoothing  as 

discussed  in  Section 

2.2 

and 

shown  in 

Figures  3-11  and 

3-12.  Because  0 is  not  a function  of  n for  the  DWD, 
frequency  resolution  is  limited  primarily  by  the  inherent 
windowing  associated  with  the  FFT.  That  is,  frequency 
resolution  is  primarily  limited  by  the  capabilities  of  the 
FFT  processor. 


3.2.2  Comparative  Examples 


Comparison  between  the  DWD  and  the  STP  can  be  best 
illustrated  through  examples.  To  establish  a like 
comparison,  equivalent  transform  lengths  and  rectangular 
windowing  will  be  applied  to  both  representations  and 
respective  magnitude  and  square  magnitude  plots  for  each 
will  be  given  in  the  figures.  The  signals  will  be  carried 
at  lower  frequencies  than  in  practical  systems.  This  is 
done  for  illustrative  clarity  and  should  not  affect  the 
purpose  of  the  examples. 
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Figures  3-13  and  3-14  give  the  respective  64  point  and 
256  point  STPs  for  a sine  pulse  on  a 400Hz  carrier  (Figure 
3-15).  For  the  sine  pulse  there  exists  no  sudden  change  in 
signal  amplitude  to  severely  degrade  the  STP  and  for  the  64 
point  transform  the  quasi-stationary  assumption  appears  to 
hold  fairly  well.  Naturally  higher  frequency  resolution  is 
attained  with  the  larger  256  point  transform,  however;  the 
quasi-stationary  assumption  over  the  window  does  not  hold 
as  well,  as  becomes  apparent  in  the  temporal  amplitude 
variation  of  the  signal.  Figures  3-16  and  3-17  are  the 
corresponding  64  and  256  point  DWDs  of  the  sine  pulse 
modulation . 

Although  the  two  examples  above  graphically  illustrate 
the  potential  of  the  DWD  for  representing  a signal  over  the 
t-o)  plane,  the  importance  of  the  DWD  has  been  its 
application  to  time-varying  signals  or  signals  whose 
dispersion  characteristics  undergo  significant  change  over 
time.  An  obvious  example  is  the  case  of  FM  signals  where 
concentration  of  the  signal  energy  shifts  over  time.  The 
larger  the  dispersion  index  the  larger  the  variance  of  the 
corresponding  spectral  output.  The  dispersion  index,  g,  is 
defined  as  the  product  of  the  period  of  the  modulating 
signal  and  its  amplitude  [35]  (i.e.  the  maximum  frequency 
deviation  from  the  carrier  frequency  of  the  FM  signal,  Af). 
For  narrowband  signals  6 should  be  on  the  order  of  <.2;  as 
0 becomes  large  (-»®)  the  resolution  band  approaches  the 
peak-to-peak  deviation  in  the  FM  signal. 


Figure  3-9.  DWD  of  square  pulse  modulation 

example  for  a 64  point  transform 
length . 


Figure  3-10. 


DWD  of  square  pulse  modulation 
example  for  a 256  point  transform 
length . 
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Figure  3-11.  DWD  of  square  pulse  modulation 

example  for  a 64  point  transform 
length  with  truncated  Gaussian 
windowing . 


Figure  3-12.  DWD  of  square  pulse  modulation 

example  for  a 256  point  transform 
length  with  truncated  Gaussian 
windowing . 
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Figure  3-13. 


STP  of  sine  pulse  modulation 
example  for  N=64. 


Figure  3-14. 


STP  of  sine  pulse  modulation 
example  for  N=256. 


51 


VALUES 


2048 
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ENO  =>  3.520  SEC 


Figure  3-15.  Sine  pulse  envelope  on  a 400Hz 
carrier . 
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Figure  3-16. 


DWD  of  sine  pulse 
example  for  N=64. 


modulation 


Figure  3-17. 


DWD  of  sine  pulse  modulation 
example  for  N=256. 
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A standard  example  of  FM  modulation  is  the  case  of 
sinusoidal  FM  on  a carrier  where  the  dispersion  index 
equates  to  that  of  the  modulation  index.  Figure  3-18 
illustrates  a sinusoidal  FM  signal  on  a carrier,  for 
Af=400Hz,  over  one-half  the  period  of  the  modulating 
signal.  The  modulation  index  is  large,  g=800,  and  the 
signal  becomes  distorted  at  the  peak  points.  For  a 
modulation  of  this  magnitude  the  spectral  band  approaches 
the  peak-to-peak  deviation.  Figure  3-19  illustrates  the 
STP  of  the  sinusoidal  FM.  Application  of  the  window  at 
each  time  slice  has  effectively  reduced  the  dispersion 
index  (notice,  the  modulation  index  and  the  dispersion 
index  are  no  longer  equivalent);  still  the  frequency  band 
at  each  time  slice  is  large.  A smaller  transform  length 
will  further  narrow  the  band;  however,  will  increase  the 
effective  bandwidth  (i.e.  reduce  resolution).  This 
translates  to  the  well-known  trade-off  associated  with  the 
STP.  That  is,  selection  of  the  window  length  with  respect 
to  frequency  resolution  and  assumption  of  stationarity . In 
the  DWD  kernel  Af  is  the  maximum  deviation,  over  the 
window,  of  the  derivative  of  [ m( t ) -m( -t ) ] , for  m(t)  the 


modulating  signal,  and  S 


DWD 


is  effectively  reduced  to 


PsTp“^ ^ over  T of  m(T-t)}  which  implies 
Figure  3-20  illustrates  the  superiority  of 
estimating  the  sinusoidal  FM. 


the  DWD  in 
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Figure  3-18.  Sinusoidal  FM  signal  for 

f =400Hz,  8=800. 
c 
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Figure  3-19. 


STP  of  sinusoidal  FM  example 
for  N=256. 
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Figure  3-20.  DWD  of  sinusoidal  FM  example 
for  N=256 . 


56 


VALUES 


Figure  3-21.  Cosine  chirp  signal. 
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Figure  3-22. 


STP  of  cosine  chirp  example 
for  N=256 . 


Figure  3-23.  DWD  of  cosine  chirp  example 
for  N=256 . 
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Figure  3-24. 


STP  of  square  pulse  FM  example 
for  N=256. 


Figure  3-25. 


DWD  of  square  pulse  FM  example 
for  N=256. 
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Figure  3-26.  Square  pulse  FM  signal  for 

f =0,  f =300Hz. 
c in 
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An  important  signal  type  to  many  applications  is  the 
cosine  chirp  (Figure  3-21)  where  the  signal  frequency 
varies  linearly  over  time.  Figures  3-22  and  3-23  are  the 
respective  256  point  STP  and  DWD  of  the  cosine  chirp.  Even 
for  the  same  transform  length  frequency  resolution  is 
poorer  for  the  STP;  this  results  from  the  postulation  that 
Af  is  small  over  the  length  of  the  window. 

Figures  3-24  and  3-25  illustrate  the  comparison 
between  the  STP  and  DWD  for  square  pulse  FM  (Figure  3-26). 
The  greatest  difference  occurs  at  the  point  where  the 
signal  spectrum  suddenly  shifts:  in  the  STP  the 
quasi-stationary  assumption  obviously  does  not  hold  and 
primary  lobes  are  present  at  both  frequencies;  in  the  DWD 
the  points  of  transition  are  represented  by  an  intermediate 
lobe  to  display  the  temporal  variation  of  the  spectrum. 

3 . 3 Smoothed-pseudo  Discrete  Wigner  Distribution 

An  important  property  of  the  WD,  and  equivalently  the 
DWD,  is  its  bilinear  nature.  That  is,  the  DWD  of  the  sum 
of  two  signals  does  not  produce  the  sum  of  the  two  DWDs  of 
the  signals  but  the  sum  of  the  DWDs  of  the  two  signals  plus 
a cross-DWD  [6].  More  explicitly,  if  y ( n ) =x ( n ) +z ( n ) and 
W (n,k)  and  W (n,k)  denote  the  respective  DWDs  of  x(n)  and 

X z 

z(n),  then  the  DWD  of  y(n)  expresses  as  below 
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N-1 


W 


,(n,k)=W  (n,k)+W  (n,k)+4Re{Z  x ( n+m) z* ( n-m) e j^nmk/Nj^ 


m=0 


(3-20) 


The  bilinear  nature  of  the  DWD  is  especially  important 
to  estimation  of  multi-component  signals  as  the  DWD 
produces  cross-term  components  at  each  DWD  slice.  For 
example,  consider  a signal  to  be  the  sum  of  two  sinusoids 
(Figure  3-27):  the  DWD  of  the  signals  produces  the  DWD  of 
each  sinusoidal  component  plus  the  cross-DWD  of  the  two 
sinusoids  which  produces  frequency  components,  called 
cross-terms,  at  dc  and  the  sum  and  difference  frequencies 
of  the  two  sinusoids.  Figures  3-28  and  3-29  are  the 
respective  STP  and  DWD  of  the  composite  800  and  600Hz  tones 
in  Figure  3-27:  the  DWD  is  corrupted  by  the  sum  and 
difference  frequencies  which  have  been  isomorphically 
mapped  down  by  a factor  of  two  due  to  oversampling  of  the 
signal.  (Note  should  be  made  that  (3-29)  is  a magnitude 
plot . ) 

The  physical  significance  of  these  cross-term 
components  is  open  to  interpretation  and  further  research. 
Chester  [36]  proposed  a significance  to  the  cross-term 
components  for  the  case  of  speech  signals  by  modelling  the 
auditory  properties  of  the  ear  as  a receiver.  For  closely 
spaced  tones  the  DWD  cross-terms  may  be  interpreted  as  a 
physical  indication  of  the  phase-coupling  between  the  two 
tones  which  would  invariably  maintain  a nonstationary 
characteristic  over  time. 
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Figure  3-27.  Sum  of  two  sinusoids  for 
fj^  = 600Hz,  f2  = 800Hz. 
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Figure  3-28.  STP  of  composite  sinusoid  example 
for  N=256 . 


Figure  3-29.  DWD , for  N=256,  of  composite 
sinusoid  example  where 
additional  cross-term 
components  at  ( £^^  + £2)/^  and 
(f^-f2)/2  are  present. 


64 


A viable  argument  as  to  the  physical  significance  of 
these  cross  terms  can  be  made  by  recalling  from  Chapter  2 
that  the  Wigner  distribution  attains  negative  values  and 
observing  that  these  negative  contributions  occur 
exclusively  in  the  cross-term  components.  To  see  this, 
consider  the  Fourier  decomposition  of  the  signal  x(t)  at 
time  t as  shown  below 

x(t)~EC  e^%^;  (3-21) 

m 
m 

from  its  bilinear  property  the  WD  of  (3-21)  equates  to  the 
sum  of  the  WDs  of  each  component  plus  the  sum  of  the 
cross-WDs  as  follows; 


W(t,w)~E|Cjjj|^5(w-Wj^)  + Ei:|Cjj^Cjjj*|cos[  (w^-“j„: 

m mn 


t]  S(  oo-(  00j^+Wj^)/2  ) . 

(3-22) 


The  first  term  in  (3-22)  is  strictly  positive  whereas  the 
second  summation  term,  corresponding  to  the  cross-term 
components,  contains  both  negative  and  positive  terms. 
Recall,  in  Section  2.2  these  negative  values  were 
demonstrated  as  important  to  the  simultaneous  preservation 
of  the  characterizing  time  and  frequency  domain  properties 
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of  the  signal.  Accordingly,  the  existence  of  the 
cross-terms  may  prove  important  to  certain  signal  types. 
For  many  applications,  however,  these  cross-terms  serve 
primarily  to  obfuscate  the  estimation  procedure.  An 
obvious  example  is  the  case  of  multipath  radio  signals 
(Sections  3.3.2  and  3.3.4). 

Flandrin  [9]  introduced  a smoothed-pseudoDWD  to  combat 
the  cross— term  problem  of  the  DWD  where  smoothing  was 
performed  on  the  kernel  as  below. 


N-1 

spW( n , k ) =2  E 1 h( m) 
m=0 


2 


g-2  j iimk/N 
m 


M-1 

E g(  m' 
' =-M+l 


) i^(m+m' ) , 


(3-23) 


where  g(m)  constitutes  smoothing  in  the  time  direction  and 
h(n)  in  the  frequency  direction.  Figure  3-30  illustrates 
the  256  point  spDWD  of  the  composite  sine  waves  pictured  in 
Figure  3-27  for  M=32.  The  value  of  M is  always  smaller 
than  N;  its  actual  value  is  system  and  application 
dependent . 

At  this  point  two  additional  points  concerning  the 
bilinearity  of  the  DWD  should  be  mentioned.  First,  the  dc 
corruption  associated  with  the  DWD  of  real  signals  may  be 
treated  as  cross-term  contribution  between  negative  and 
positive  frequency  components.  Application  of  the  analytic 
signal  can  then  serve  to  eliminate  the  dc  corruption. 
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Figure  3-30.  Smoothed-psuedo-DWD  of  composite 
sinusoid  example  for  M=32  and 
N=256 . 
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Second,  from  Equation  (3-20)  it  follows  the  DWD  is  easily 
degraded  when  noise  is  present.  That  is,  if  z(n) 
represents  additive  noise  then  at  each  each  DWD  slice,  n, 
the  noise  contribution  comprises,  in  addition  to  the  DWD  of 
the  noise,  a cross-DWD  of  the  signal  with  the  noise  defined 
by  the  DFT  of  x^(  m)  z^*  ( -m) +z^^(  m)  x^*  ( -m)  . For  this  reason 
the  DWD  has  been  treated  as  an  estimator  that  requires 
clean  signals  [36].  This  constraint  severely  limits  the 
DWDs  application  potential.  Fortunately  the  spDWD  can  be 
used  to  effect  significant  noise  reduction. 

The  spDWD  can  be  interpreted  as  the  DSP  analog  to  the 
weighted-WD  of  Section  2.2.2,  where  smoothing  is  effected 
on  each  pDWD  slice  as  below 


M-1 

spW(n,k)=Z  g(m)W(n+m,k) . 
m=-M+l 


(3-24) 

the  spDWD 

can 

be 

of  the  DWD 

where 

the 

preserve 

the 

moment  properties  of  the  DWD  the  weighting  functions  are 
real  symmetric.  From  this  interpretation  follows  that  the 
spDWD  weighting  functions  can  be  interpreted  as  filtering 
and  modulation  operations  on  the  DWD  as  follows  (Section 
2.2.2  and  2.2.3) 
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spW( n , k ) 


M-1  N-1 

= E E W(m,i)H' (k-i)g(n-m) ; 
m=-M+l  i=0 


H' (k)=E^h^(n)e  j^nkn/N^  (3-25) 
n=0 

In  this  form  the  wealth  of  knowledge  on  filtering  and 
modulation  may  be  applied  to  the  spDWD  to  achieve  superior 
estimation  for  the  two  primary  conditions  under  which  the 
DWD  produces  poor  results;  multicomponent  signals  and  noisy 
environments.  Two  primary  advantages  can  be  seen  from  the 
separability  of  the  spDWD  weighting  function 
e(n,k)=g(n)HMk) . First,  the  avoidance  of  time/frequency 
smoothing  trade-off  associated  with  the  STP  and  second,  the 
ability  to  effect  smoothing  similar  to  that  performed  on 
the  STP  (to  improve  performance)  without  the  computational 
burden  of  computing  multiple  FFTs. 

3.3.1  Examples  of  Cross-term  Reduction 

Consider  for  example  the  sum  of  two  closely  spaced 
cosine  chirps,  as  pictured  in  Figure  3-31.  The  DWD  (Figure 
3-32)  produces  sufficient  resolution  between  the  two  chirps 
but  generates  an  additional  cross-term  component  between 
the  lobes  at  each  time  slice.  Application  of  the  spDWD 
(Figure  3-33)  temporally  smooths  the  signal  which 
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effectively  filters  out  the  cross-term  component  as  the 
distance  between  the  signal  lobes  increases.  For  the  case 
where  the  signal  lobes  are  closely  spaced  a larger  value  of 
M is  needed  to  effect  a narrow  band.  However  for 
time-varying  signals,  M must  be  small  enough  such  that  the 
convolution  over  which  the  signals  inner  product  and  the 
temporal  weights  are  taken  does  not  exceed  the  local  time 
of  s tationar i ty ; otherwise  the  nonstationary 
characteristics  of  the  signal  become  obscured  from  too  much 
temporal  smoothing.  Figure  3-34  is  the  spDWD  of  the  dual 
cosine  chirp  for  the  case  of  smoothing  over  the  entire 
window;  although  the  cross-terms  have  been  eliminated, 
temporal  smoothing  has  degraded  the  signals  time-varying 
characteristic  and  the  spDWD  representation  approaches  that 
of  the  STP  (Figure  3-35)  which  postulates  stationarity  over 


the  window. 

The  effect  on  the  time-varying  characteristic  of  a 
signal  relative  to  the  choice  of  M becomes  apparent  when 
the  spDWD  is  applied  to  the  pulse  FM  signal  of  Figure  3-21 
for  increasing  values  of  M (Figures  3-36,  3-37,  and  3-38). 
For  M small  (=4)  the  cross-term  components  are  effectively 
eliminated.  Recall,  the  dc  corruption  can  be  treated  as 
cross-term  components  between  negative  and  positive 
frequencies;  therefore,  if  the  spectral  lobes  are 
sufficiently  away  from  dc,  minimal  smoothing  will  eliminate 
the  dc  terms.  As  M is  increased  the  temporal 
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Figure  3-31.  Dual  cosine  chirp. 
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Figure  3-32.  DWD  of  dual  cosine  chirp  example 
for  N=256. 


Figure  3-33. 


Smoothed-psuedo-DWD  of  dual 
cosine  chirp  for  M=12  and 
N=256 . 
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Figure  3-34.  Smoothed-psuedo-DWD  of  dual 
cosine  chirp  for  M=128  and 
N=256  . 
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Figure  3-35. 


STP  of  dual  cosine  chirp. 
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Figure  3-36.  Smoothed-psuedo-DWD  of  square 
pulse  FM  example  for  M=4  and 
N=256 . 


Figure  3-37.  Smoothed-psuedo-DWD  of  square 
pulse  FM  example  for  M=32  and 
N=2 56  . 
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Figure  3-38. 


Smoothed— psuedo-DWD  of  square 
pulse  FM  example  for  M=128  and 
N=256 . 
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characteristic  of  the  spDWD  approaches  that  of  the  STP, 
especially  at  the  transient  points  where  the  local  time  of 
stationarity  approaches  zero  (notice  at  these  points  any 
amount  of  smoothing  exceeds  the  local  time  of  stationarity 
for  the  pulse  FM ) . 

3.3.2  Multipath  Radio  Example 

Figure  3-39  is  a cosine  chirp  with  an  attenuated 
reflection  delayed  in  time  and  shifted  in  frequency.  A 
signal  of  this  type  is  called  a multipath  signal  and 
presents  a common  problem  in  radio  transmission  between  low 
flying  aircraft.  The  attenuated  component  results  from 
ground  reflection  and  interferes  with  the  direct  path 
signal.  If  the  reflection  contribution  can  be  properly 
detected  its  presence  may  be  accounted  for  at  the  receiver. 

Figures  3-40  and  3-41  give  the  respective  256  point 
STP  and  DWD  of  the  multipath  signal  in  Figure  3-39.  In  the 
STP  existence  of  the  echo  cannot  be  distinguished  whereas 
in  the  DWD  three  separate  tones  appear:  the  tone 
corresponding  to  the  original  chirp  the  tone  corresponding 
to  the  echo  whose  DWD  is  shifted  in  both  frequency  and  time 
from  the  DWD  of  the  original  chirp  and  a cross-term  between 
them.  Application  of  the  spDWD  (Figure  3-42)  can  be  used 
to  eliminate  that  part  of  the  cross-term  that  results  from 
the  two  frequency  components;  however,  the  contribution  to 
the  cross-term  component  that  results  from  the  delay 
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Figure  3-39.  Multipath  radio  example:  cosine 
chirp  envelope  with  attenuated 
reflection  delayed  in  time  and 
shifted  in  frequency  where 
t^=T/6  and  f^=200Hz. 
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Figure  3-40,  STP  of  multipath  chirp  example 
for  N=256. 


Figure  3-41.  DWD  of  multipath  chirp  example 
for  N=256. 
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Figure  3-42. 


Smoothed-psuedo-DWD  of  multipath 
chirp  example  for  M=32  and  N=256. 
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between  the  chirps  is  relatively  invariant  to  the  smoothing 
procedure.  That  is,  recall  from  Equation  (3-25)  the 
smoothing  procedure  effects  a filtering  operation  on  the 
signal  at  each  time  slice,  or  a moving  filter.  In 
addition,  in  the  smoothing  procedure,  the  shape  of  the 
delayed  pulse  becomes  distorted  as  the  spDWD  performs 
temporal  smoothing  on  the  signal  characteristics. 

Consider  for  example  the  sum  of  a cosine  chirp  and  a 
delayed  version  of  the  chirp  with  no  frequency  shift 
(Figure  3-43).  In  the  DWD  (Figure  3-44)  three  tones  are 
present.  The  center  lobe  is  a cross-term  produced  by  the 
two  outer  lobes  which  result  from  the  time  displacement  in 
the  DWD  of  each  chirp.  Care  must  be  taken  not  to  confuse 
the  time  displaced  lobes  with  frequency  shifted  lobes.  The 
difference  for  the  linear  FM  example  is  clear:  in  this 
example  the  displacement  does  not  linearly  increase  over 
time  as  in  the  DWD  of  the  dual  chirp  of  Figure  3-32. 
Figure  3-45  is  the  spDWD  for  M=12.  While  smoothing 
eliminates  the  dc  components  that  result  from  cross-terms 
between  the  negative  and  positive  frequencies  of  the  cosine 
chirp  (Figure  3-33),  the  cross-term  component  between  the 
time  displaced  chirps  cannot  be  eliminated.  Additional 
smoothing  (Figure  3-46)  does  not  reduce  the  relative 
amplitude  of  the  cross-term  but  obscures  the  output  by 
smoothing  of  the  temporal  signal  properties.  Figure  3-47 
is  the  STP  of  the  delayed  chirp  which  precludes  the 
existence  of  a delayed  chirp  altogether. 
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Figure  3-43.  Chirp  plus  delay  where  x^=T/6 . 
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Figure  3-44. 


DWD  of  chirp  plus  delay  example 
for  N=256. 


Figure  3-45.  Smoothed-psuedo-DWD  of  chirp 
plus  delay  example  for  M=12 
and  N=256 . 
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Figure  3-46.  Smoothed-psuedo-DWD  of  chirp 
plus  delay  example  for  M=64 
and  N=256. 


Figure  3-47.  STP  of  chirp  plus  delay  example 
for  N=256 . 
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Recall  the  double  Fourier  transform  relationship 
between  the  DWD  and  the  ambiguity  function  discussed  in 
Chapter  2.  In  the  ambiguity  function  instead  of  the  signal 
waveform,  the  delay  and  doppler  shift  of  the  signal  are 
ideally  represented  by  location  of  a secondary  lobe.  (The 
classes  of  practical  ambiguity  diagrams  and  their 
limitations  are  well  documented  [33]  and  their  inclusion  is 
not  necessary  to  this  discussion.)  A well-known  property  of 
the  AF  is  that  to  obtain  maximum  resolution  of  the 
frequency  and  time  delay  while  keeping  the  ambiguity  which 
results  from  additional  peaks  to  a minimum  requires  the 
bandwidth-time  (BT)  product  of  the  signal  be  as  large  as 
possible.  That  is,  the  worst  waveform,  i.e.  that  which 
produces  maximum  variance  in  the  AF  plane,  is  the  Gaussian. 
Following  with  the  Fourier  transform  relationship  between 
the  AF  and  the  DWD,  the  optimum  window  to  minimize  the 
variance  in  the  DWD  plane  is  a Gaussian.  By  selecting  the 
filtering  operation  to  also  be  of  Gaussian  form  (recall  the 
transform  of  the  Gaussian  produces  a Gaussian)  smoothing  of 
the  time  properties  can  be  minimized.  (A  formal  discussion 
of  Gaussian  smoothing  of  the  WD  is  provided  in  Section 
2.2.3.  ) 


Figures  3-48  and  3-49  are  the  respective  spDWDs  for 
h(n)  a Gaussian,  g(n)  rectangular,  and  for  both  h(n)  and 
g(n)  Gaussian.  A subtle  improvement  in  Figure  3-49  can  be 
observed.  Additional  improvement  can  be  effected  by 
increasing  M and  reducing  the  variance  of  g(n),  or 
effectively  increasing  the  variance  of  the  Gaussian  filter. 
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Figure  3-48.  Smoothed-psuedo-DWD  of  multipath 
chirp  example  for  M=32  and  N=256 
where  Gaussian  windowing  has  been 
applied . 


Smoothed-psuedo-DWD  of  multipath 
chirp  example  for  M=32  and  N=256 
where  Gaussian  windowing  and 
filtering  have  been  applied. 


Figure  3-49. 
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3.3.3  Examples  of  Noise  Reduction 


Figures  3-50  and  3-51  are  the  respective  STP  and  DWD 
of  a weak  tone  in  the  presence  of  additive  pseudo-random 
noise  for  a signal-to-noise  ratio  of  -lOdb  (Figure  3-52). 
In  the  STP  the  noise  greatly  obscures  the  tone.  In  the  DWD 
amplification  of  the  noise  precludes  observation  of  the 
tone.  Figure  3-53  is  the  corresponding  spDWD  for  M=32 
(g(n)  rectangular)  and  Figure  3-54  is  the  spDWD  for 
smoothing  over  the  length  of  the  window. 

Figures  3-55  and  3-56  are  the  STP  and  DWD  of  the  chirp 
plus  noise  in  Figure  3-57.  Figure  3-58  is  the  spDWD  for 
both  smoothing  parameters  weighted  as  a truncated  Gaussian 
where  truncation  occurs  at  the  variance.  By  truncating  the 
window  at  the  variance  the  signal  lobe  remains  narrow; 
otherwise,  the  signal  and  noise  resolution  would  degrade 
unless  a longer  tranform  length  was  applied.  (Recall, 
however,  the  quasi-stationary  approximation  which  limits 
the  transform  length  for  the  STP,  even  for  the  case  of 
windowing,  does  not  apply  to  the  DWD  and  a longer  transform 
length  may  be  a viable  option.)  Alternatively,  Gaussian 
weighting  of  the  filter  would  require  a larger  value  of  M 
to  sufficiently  smooth  the  noise.  However,  a larger  value 
of  M will  widen  the  lobe  as  the  temporal  characteristics 
become  smoothed  unless  the  variance  is  made  sufficiently 
small.  In  addition,  if  the  value  of  M is  limited  by  the 
spDWD  processor  truncation  can  be  effectively  employed. 
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Figure  3-50.  STP  of  sinusoid  with  additive  noise 
example  for  N=256. 


Figure  3-51.  DWD  of  sinusoid  with  additive  noise 
example  for  N=256. 
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Figure  3-52.  Sinusoid  with  additive  psuedo-random 
noise  for  a signal-to-noise  ratio 
of  -lOdb. 
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Figure  3-53.  Smoothe(d-psuedo-DWD  of  sinusoi(d 
with  a(d(ditive  noise  example  for 
M=32  and  N=256. 


Figure  3-54.  Smoothed-psuedo-DWD  of  sinusoid 
with  additive  noise  example  for 
M=128  and  N=256. 
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Figure  3-55.  STP  of  chirp  plus  noise  example 
for  N=256. 


Figure  3-56.  DWD  of  chirp  plus  noise  example 
for  N=256. 


VALUES 


Figure  3-57.  Cosine  chirp  with  additive 
psuedo-random  noise  for  a 
signal-to-noise  ratio  of  -3db. 
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Figure  3-58.  Smoothed-psuedo-DWD  of  chirp 
plus  noise  example  for  M=12 
and  N=256  where  Gaussian 
windowing  and  filtering  have 
been  applied. 
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3.3.4  Multipath  Radio  Signal  with  Receiver  Noise 


As  a final  example,  consider  the  multipath  signal  of 
Figure  3-39.  In  practical  applications  the  radio  receiver 
introduces  an  additive  zero  mean  Gaussian  noise  term  [37]. 
Figure  3-59  is  the  resulting  DWD  when  severe  receiver  noise 
(Figure  3-60)  is  added  to  the  composite  chirp  plus 
reflected  chirp  of  Figure  3-39  to  produce  the  corrupted 
signal  in  Figure  3-61.  Figures  3-62  and  3-63  are  the 
respective  spDWD  and  spDWD  with  Gaussian  smoothing  for 
M=32.  In  the  spDWD  the  presence  of  the  reflected  signal 
remains  obscured;  Gaussian  smoothing  offers  some 
improvement.  A large  value  of  M effects  too  much  smoothing 
of  the  temporal  properties  of  the  signal  and  serves  only  to 
degrade  the  resultant  output  (Figures  3-64).  However,  by 
applying  a large  value  of  M with  Gaussian  smoothing  where 
the  standard  deviation  (STD)  of  the  Gaussian  filter  is 
relatively  large  or,  equivalently,  the  STD  of  the  temporal 
weight  is  relatively  small,  satisfactory  noise  reduction 
can  be  achieved  and  the  reflection  component  becomes 
visible  (Figure  3-65).  Figure  3-66  is  the  spDWD  for  no 
Gaussian  smoothing  where  the  value  of  M is  equal  to  the  STD 
of  the  temporal  weight  in  Figure  3-65.  Here  the  reflection 
component  remains  obscured  and  both  the  direct  and 
reflected  signal  components  are  distorted. 

The  spDWD  is  sensitive  to  selection  of  the  filter  STD. 
Reducing  the  STD  of  the  Gaussian  filter  slightly  (Figure 
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3-67)  reduces  the  cross-term  components  which  have  been 
amplified  by  the  noise.  However,  selecting  too  small  a STD 
for  the  Gaussian  filter  will  produce  too  much  temporal 
smoothing  on  the  signal  as  shown  in  Figure  3-68.  Figures 
3-69  and  3-70,  the  STP  for  rectangular  and  Hamming  windows, 
have  been  added  for  the  sake  of  comparative  completeness. 
Application  of  the  Hamming  window  offers  some  improvement. 
Still,  in  both  figures  the  noise  precludes  observation  of 
the  reflection  component. 

Chapter  4 will  concentrate  on  practical  means  of 
implementing  the  WD,  the  DWD,  and  the  spDWD  for  both  real 
and  complex  signals. 
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Figure  3-59.  DWD  of  multipath  chirp  with 
receiver  noise  example  for 
N=256 . 
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Figure  3-60.  Zero  mean  Gaussian  noise  with  a 
standard  deviation  of  700. 
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Figure  3-61.  Multipath  radio  with  receiver 
noise  example:  cosine  chirp 
envelope  with  attenuated 
reflection  delayed  in  time  and 
shifted  in  frequency  where 
Tj=T/6  and  fj=200Hz  in  the 
presence  of  additive  Gaussian  noise. 
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Figure  3-62.  Smoothed-psuedo-DWD  of  multipath 
chirp  with  receiver  noise  example 
for  M=32  and  N=256. 


Figure  3-63.  Smoothed-psuedo-DWD  of  multipath 
chirp  with  receiver  noise  example 
for  M=32  and  N=256  where  Gaussian 
windowing  and  filtering  have 
been  applied. 
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Figure  3-64.  Smoothed-psuedo-DWD  of  multipath 
chirp  with  receiver  noise  example 
for  M=128  and  N=256. 


Figure  3-65.  Smoothed-psuedo-DWD  of  multipath 
chirp  with  receiver  noise  example 
for  M=128  and  N=256  where  Gaussian 
windowing  and  filtering  have 
been  applied  and  the  STD  of  the 
filter  is  18Hz . 
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Figure  3-66.  Smoothed-psuedo-DWD  of  multipath 
chirp  with  receiver  noise  example 
for  M=64  and  N=256. 
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Figure  3-67.  Smoothed-psuedo-DWD  of  multipath 
chirp  with  receiver  noise  example 
for  M=128  and  N=256  where  Gaussian 
windowing  and  filtering  have 
been  applied  and  the  STD  of 
the  filter  has  been  reduced 
to  llHz. 
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Figure  3-68. 


Smoothed-psuedo-DWD  of  multipath 
chirp  with  receiver  noise  example 
for  M=128  and  N=256  where  Gaussian 
windowing  and  filtering  have 
been  applied  and  the  STD  of 
the  filter  has  been  reduced 
to  4Hz. 
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Figure  3-69.  STP  of  multipath  chirp  with 
receiver  noise  example  for 
N=256. 


Figure  3-70.  STP  of  multipath  chirp  with 
receiver  noise  example  for 
N=256  where  Hamming  windowing 
has  been  applied. 


CHAPTER  4 


WIGNER  IMPLEMENTATIONS 


4 . 1 Wigner  Processing  Using  Standard  FFTs 


4.1.1  General  Wigner  Processing 


An  attractive  feature  of  the  DWD  is  that  it  can  be 
computed  via  the  DFT  of  an  inner  product  function  (or  for 
the  spDWD  a weighted  inner  product  function)  defined  by  the 
the  inner  product  of  the  signal  shifted  in  time  and  its 
complex  conjugate  reversed  in  time  and  referred  to  as  the 
DWD  kernel.  In  this  form  the  DWD  (and,  equivalently  the 
spDWD)  readily  lends  itself  to  FFT  implementations.  In 
short,  a DWD  processor  comprises  two  major  subsystems: 
kernel  generation  and  a DFT. 

Kernel  generation  for  the  DWD  is  relatively 
straightforward;  windowing  is  performed  prior  to 
computation  of  the  inner  product.  The  input  sequence  is 
loaded  into  an  N point  buffer.  Once  the  buffer  is  loaded. 
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one  half  of  the  DWD  kernel  is  generated  by  multiplying  the 
initial  data  point  of  the  input  sequence  x(0)  with  the  last 
point  of  the  input  sequence  x(N-l)  to  form  the  first  point 
of  the  DWD  kernel,  multiplying  x(l)  with  x(N-2)  to  form  the 
second  point  etc.,  until  the  kernel  from  0 to  N/2-1  is 
completed.  For  real  data  the  DWD  kernel  is  symmetric  about 
its  center  and  each  product  forms  two  kernel  values.  For 
complex  data  the  kernel  is  congugate  symmetric  about  the 
center.  Durrani  et  al.  [38]  presented  a systolic 
architecture  for  small  N DWDs  of  real  data.  As  discussed 
in  Chapter  2,  for  real  input  data  rearrangement  is  not 
actually  needed  in  the  DWD  computation. 

In  Section  3.3  the  spDWD  was  advanced  and  shown  to 
maintain  several  attractive  features;  especially  for  the 
two  primary  conditions  under  which  the  DWD  produces  poor 
results;  multicomponent  signals  and  noisy  environments. 
Despite  the  advantages  of  the  spDWD,  the  increased 
computational  requirements  associated  with  producing  the 
spDWD  kernel  are  clear  from  its  definition.  Flandrin  and 
Martin  proposed  a straightforward  hardware  implementation 
for  the  spDWD  capapable  of  "quasi-real-time"  data  rates 
[39]  . 

This  section  will  present  a high  speed  concurrent 
architecture  composed  of  single  task  processor  cells  to 
generate  the  spDWD  kernel.  Each  cell  operates 
independently;  once  all  inputs  (or  operands)  have  been 
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supplied,  the  processor  task  is  performed  then  output  as  an 
operand  to  the  other  adjacent  cell.  For  illustrative 
purposes,  an  N=4,  M=3  point  processor  is  given  in  Figure 
4-1.  The  first  set  of  cells  perform  complex  conjugate 
multiplication  on  the  input.  The  output  products,  i(m,n), 
from  these  cells  constitute  as  operands  to  the  complex 
multiply-add  cells.  For  each  multiply-add  cell,  the 
operand  i(m,n)  is  shifted  in  and  multiplied  by  the  temporal 
weight  g(m)  which  has  also  entered  the  cell.  The  product 
is  then  added  to  the  composite  output,  or  sum  of  products 
elements  that  have  entered  the  cell  from  the  adjacent 
^^^biplier  cell  to  form  the  new  sum  of  products  element, 

^ be  input  to  the  other  adjacent 
cell.  After  an  intial  startup  delay  of  2M+1  complex 
^^Itiply  operations,  the  total  smoothed  WD  inner  product 
output  for  each  consecutive  time  slice  is  supplied  to  the 
DFT  buffers  after  each  operation  period.  The  startup  delay 
is  minimal  compared  to  the  savings  afforded  by  having 
continuous  output  after  the  initial  startup. 

In  most  applications  M is  typically  much  smaller  than 
N to  avoid  masking  the  nonstationary  characteristics  of  the 
signal  and  longer  length  spDWDs  are  computed  by  setting 
N=N'N''  and  completing  several  passes  through  the  processor 
to  effect  the  total  WD  smoothed  inner  product  output.  (By 
removing  the  weighted  terms  at  m=0  and  feeding  the 
composite  product  into  S^,  temporal  smoothing  may  be 
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g(-2)  g(-l) 


8(0)  g(l)  g(2) 


Figure  4-1. 


N=4,  M=3  smoothed-psuedo-DWD 
processor . 
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similarily  decomposed.  Each  weighted  term  at  m=0  is  then 
added  prior  to  windowing.)  For  real  input  data  the 
processing  task  comprises  all  real  multiplies  and  is 
greatly  simplified.  As  with  the  DWD  kernel,  the  weighted 
kernel  is  symmetric  for  real  data  and  conjugate  symmetric 
for  complex  data  about  its  center. 

Naturally,  the  selection  of  FFT  for  the  DFT  stage  of 
of  the  Wigner  processor  is  both  system  and  application 
dependent.  Reported  work  on  the  FFT  is  voluminous  and  an 
extensive  reference  listing  of  efficient  FFT  algorithms  is 
compiled  in  [40].  Detailed  tutorials  on  radix-r  type  FFT 
decompositions  can  be  found  in  [41,42]  and  in  more  recent 
texts  which  concentrate  on  number  theoretic  transforms 
[43,44].  A collection  of  fundamental  papers  on  NTTs  is 
included  in  [43].  A firm  foundation  of  the  modern  FFT  is 
provided  by  Winograd  [45-47].  Efficient  FFT  FORTRAN 
programs  are  provided  in  [48]. 

Power-of-two  algorithms  generally  require  more 
computations  than  the  number  theoretic  transforms;  however 
the  regularity  of  the  butterfly  structure  and  the  ability 
to  perform  in  place  computations  has  made  the  radix-2  and 
radix-4  algorithms  a popular  choice.  Recently  the 
split-radix  FFT  has  become  popular  because  it  reduces  the 
number  of  additions  and  nontrivial  multiplications  while 
maintaining  the  regularity,  flexibility,  and  ability  to 
perform  in  place  computations  associated  with  the 
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power-of-two  algorithms.  The  split 
decomposition  such  that  a radix-2  map 
even  indexed  terms  and  a radix-4  map 
indexed  terms.  The  split-radix  is  det 
efficient  programs  are  provided. 


radix  FFT  is  based  on 
is  applied  to  the 
is  applied  to  the  odd 
ailed  in  [49-51]  and 


4.1.2  Computational  Reduction  for  Real  Data 

The  major  computational  burden  associated  with 
producing  either  the  DWD  or  spDWD  of  real  data  sequences 
has  traditionally  been  the  FFT  stage  of  processing. 
Chester  [36]  proposed  a DWD  analyzer  for  which  the  inner 
product  and  associated  windowing,  data  rearrangement,  etc. 
was  computed  via  software  then  applied  to  an  FFT  that  was 
constructed  using  special  purpose  hardware.  In  his 
implementation  Chester  used  the  symmetry  property  of  the 
DWD  inner  product  to  compute  two  DWDs  during  one  FFT  cycle. 
Consequently,  Chester  concluded  that  the  DWD  has  nearly 
twice  the  throughput  of  the  STFT  because  two  DWDs  can  be 
computed  via  one  FFT.  A counter  to  this  claim  has  been 
that  for  real  data  an  N point  DFT  can  be  computed  via  an 
N/2  point  FFT  because  the  corresponding  output  is  conjugate 
symmetric  [41].  Using  the  symmetry  of  both  the  DWD  inner 
product  input  and  the  corresponding  DWD  output,  two  DWDs 
can,  in  fact,  be  computed  via  one  N/2  point  FFT. 

To  see  this  let  ij^(m)  and  i2(ni)  each  denote  respective 
At  the  output  stage  of  the  kernel  generator 


DWD  kernels. 
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each  kernel  sequence  is  redefined  as  an  N/2  point  sequence 
as  follows: 


r^(m)  = rj^g(m)  + rj^^(m)  ;m=0,l,  . . . , N/2-1 ; m=l , 2 , (4-la) 


where 


;n=l,2 


(4-lb) 


are  even  funtions  and 


^no^  2m+l ) -ij^(  2m-l ) ;n=l , 2 


(4-lc) 


are  odd  functions.  The  DFT  of  n=l,2  expresses  as 

follows : 


R (k)=R  (k)+R  (k);n=l,2 

n ne  no 


(4-2) 
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where  R„^(k)  and  (k)  are  the  respective  DFTs  of  r (m) 
ne  no  ne 

and  r^^(m)  for  n=l,2.  The  first  term  R^e  ^ ^ equivalent 

to  the  DFT  over  m even  of  the  original  sequence 

below: 


N/2-1 
= 2/N  E i 
ra=0 


n 


N-1 

=2/N  E i (m)e 
m=0  ^ 
m even 


- j 2 iimk/N 


n=l , 2 . 


(4-3) 


The  second  term  reduces  as  follows: 


Rno(M  = 2/N 


N/2-1 
2 [i 
m=0 


n 


(2m+l)-i^(2m-l)  ]e 


=2/NE  i (m)e  32n(m  l)k/N 
m=0 
m odd 

-2/NU^i  (m)e"^^'^^"'‘^^  ^’^/'^;n=l,2. 
m=0 
m odd 


-2/N(  e-3  2 j 2 Mk/N" j,/  I I j 2 nmk/N  j , 

m=0  ^ (4-4) 

m odd 
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Now,  by  applying  r^(m)  to  the  real  input  of  the  N/2 
point  FFT  and  r2(ni)  to  the  imaginary  part  of  the  N/2  point 
FFT,  the  effective  input  is  the  following: 


y(n)=r^(m)+jr2(m);m=0,l,... ,N/2-l 


(4-5) 


with  corresponding  output, 


Y(  k )=Rie(  k)  + jRio(  k ) + jR2e(  1<  )+R2o^  ) ; 


k=0 , 1 , . . . ,N/2-l , 


(4-6) 


by  definition  of  the  DFT.  Applying  the  folding  properties 
at  the  output  of  the  N/2  point  FFT  (4-6),  the  DWDs  W(n,k) 
n=l,2  of  the  original  kernel  sequences  ifj(ni)  n=l,2  can  be 
obtained  from  the  sum  of  the  DFT  over  m even  and  the  DFT 
over  m odd  which  via  Equations  (4-3)  and  (4-4)  equate  to 
the  following: 


W( 1 , 0 )=l/4Re{ Y( 0 )+Y(N) }+l/4lm{ Y( 0 )-Y( N) } , 


W( 2 , 0 )=l/4lm{ Y( 0 )+Y(N) }+l/4Re {-Y( 0 )+Y(N) } ; 


(4-7a) 
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for  k=l ,2 , . . . , N/4 , 


W(  1 ,k  )=l/4Re{Y(  k )+Y(N-k  ) } +1/4 Im^ ^ Y(N  k)  } 

2sin ( 2 nk/N ) ' 


W(2,k)=l/4lm{Y(k)+Y(N-k)  } +l/4Re  ^ ^ ^ ^ ^ . 

2sin( 2nk/N ) ' 

(4-7b) 


for  k=N/4+l , . . . , N/2-1 , 

W(l,k)=l/2Re{Y(N/2-k) } , 

W(2,k)=l/2lm{Y(N/2-k) } ; (4-7c) 

and , 

W( 1 ,N/2 )=l/4Re{ Y( 0 )+Y(N) } +1/4 Im{ Y ( 0 ) -Y( N ) } , 

(4-7d) 

W(2,N/2)=l/4lm{Y(0)+Y(N) } +1/4 Im{ -Y ( 0 ) +Y ( N ) } ; 
where  for  real  input  the  DWD  is  symmetric  about  N/2 . 


Ill 


4 . 2 High-speed  Processing  Via  Algebraic  Mappings 


4.2.1  Problem  Statement 


For  many  applications  restriction  to  real  signals  is 
not  practical.  In  communication  systems,  for  example, 
demodulation  techniques  at  the  receiver  produce  complex 
encoded  data.  In  addition,  application  of  an  analytic 
signal  can  be  shown  to  have  certain  advantages  specifically 
for  DWD  processing,  several  of  which  are  listed  below; 

elimination  of  information  redundancy  associated  with 
real  sequences; 

elimination  of  cross  terms  that  result  from 
interference  between  positive  and  negative  frequencies 
and  degrade  dc  information; 

the  effective  bandwidth-time  product  is  reduced; 

sampling  at  the  Nyquist  rate  is  sufficient  to  prevent 
frequency  aliasing  of  the  DWD  (for  real  signals 
sampling  must  be  at  least  double  the  Nyquist  rate); 
and , 

instantaineous  frequency  can  be  obtained  from  the  DWD. 

Whether  real  or  complex  data  are  to  be  used,  computation  of 
the  DWD  is  explicitly  complex  multiply  intensive  and  for 
the  DWD  to  be  a viable  instrument  in  real  time  applications 
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requires  a processor  capable  of  performing  complex 
operations  at  very  high  data  rates. 

In  general,  to  complete  a complex  multiply  requires 
two  real  multiplies  and  two  adds  and  traditional  methods  of 
performing  complex  arithmetic  operations  are  both  hardware 
complex  and  of  long  latency.  To  overcome  this  problem 
alternative  number  systems  have  been  explored.  The  residue 
number  system  (RNS),  due  to  its  potential  for  high-speed 
computations,  has  received  much  attention.  In  the  residue 
number  system  each  integer  is  represented  via  a set  of 
smaller  wordlength  integers,  called  residues. 

The  attractive  feature  of  the  RNS  is  that 
multiplication  (and/or  addition)  can  be  completed  by 
independently  performing  multiplication  (and/or  addition) 
on  each  of  the  smaller  wordlength  residues.  Therefore,  in 
the  RNS  each  operation  can  be  performed  via  table-lool<up 
using  high-speed  semiconductor  memory.  Numerous  papers 
have  been  published  on  application  of  the  RNS  to  transform 
design  [52-57].  In  most  cases,  however,  the  RNS  has  not 
achieved  a distinct  advantage  over  conventional 
architectures.  The  principle  obstacle  to  achieving  an 
advantage  has  been  a high  scaling  and  residue  to  decimal 
conversion  overhead,  whether  using  multi-moduli  Chinese 
remainder  theorem  (CRT)  or  mixed-radix  conversion  ( MRC ) 
routines.  That  is,  the  RNS  is  an  integer  system;  unless 
some  form  of  magnitude  scaling  is  performed  at  intermediate 
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stages  of  the  processing  task  the  dynamic  range  extension 
requirements  become  inordinate.  Unfortunately,  magnitude 
scaling  requires  conversion  from  RNS  which  is  a cumbersome 
operation.  This  fact  has  limited  the  successful 
application  of  the  RNS  in  DSP. 

If  the  moduli  used  to  produce  each  residue  satisfy 
certain  constraints  a modular  (residue)  system  called  the 
quadratic  residue  numbering  system  (QRNS)  is  defined.  In 
the  QRNS  an  isomorphism  exists  for  which  the  number  of  real 
multiplies  and  adds  necessary  to  perform  a complex  multiply 
can  be  reduced  by  two  and  one  respectively.  For  this 
reason,  despite  the  scaling  overhead  associated  with  RNS 
based  implementations,  application  of  the  QRNS  to  complex 
multiply  operations  is  warranted.  Recently  Taylor  [58] 
introduced  a single-modulus  complex  arithmetic  logic  unit 
that  utilized  the  QRNS.  Taylor's  SM-QRNS  reduced  the 
previous  obstacle  of  multi-moduli  RNS,  namely 
residue-to-decimal  conversion,  down  to  a set  of  trivial 
shift-adds  [58].  As  a result,  both  fast  and  hardware 
elegant  computational  SM-QRNS  primitives  could  be  designed. 

A smoothed-pseudo-DWD  processor  implementation  was 
shown  to  be  well  suited  to  the  SM-QRNS  [59].  The  DWD  was 
shown  to  maintain  certain  properties  that  make  it 
especially  attractive  to  the  QRNS  and  a reduction  in 
computation  time  for  the  SM-QRNS  was  mathematically 
conjectured.  The  question  of  nested  multiplies  in  the  QRNS 
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environment  end  essociated  magnitude  scaling  overhead  was 
advanced  and  alernative  FFT  structures  to  reduce  nested 
multiplies  were  proposed  for  the  DFT  stage  of  the  DWD 
processor . 

In  this  section,  a high  performance  SM-QRNS  DWD 
estimator  that  uses  standard  FFT  structures  will  be  shown 
to  be  a practical  reality.  DWD  estimators  for  both  radix-2 
and  radix-4  DIT  FFT  subsystems  will  be  provided  and 
quantified  in  terms  of  increase  in  total  system  throughput 
relative  to  transform  length  and  wordsize.  The  SM-QRNS 
will  be  shown  to  reduce  the  effective  time  delay  of  the 
spDWD  kernel  generator  proposed  in  Section  4.1  for  the  case 
of  complex  data.  Because  the  temporal  smoothing  parameter 
is  optimally  real-symmetric  the  SM-QRNS  kernel  generation 
in  [59]  will  be  modified  to  reduce  the  hardware  for  longer 
wordlengths.  First,  the  principles  of  the  QRNS  and  its 
corresponding  isomorphism  for  multiply  reduction  will  be 
developed . 

4 • 2 . 2 Quadratic  Residue  Number  System  Principles 


Given  a finite  set  of  integers,  each  integer  asZ 


M 


has  a unique  RNS  representation  given  by: 


a->(  a^^ , a2  , . . . , a^  ) ; a^=amodp^eS 


k^^k 


where 
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M=n  p,  ; gcd(p.  ,p  . )=1. 
k = l ^ ^ 


(4-8) 


Signed  numbers,  ae [ -M/2 , M/2 ) , are  represented  as  in  (4-8) 

if  a>0;  if  a<0,  then  aj^=(M- 1 a | )modpj^ . The  residue  sets 

( Zj^ , + , ' , 0 , 1 ) , k=l,2,...,L  each  form  a field  and  integer 

addition  and  multiplication  can  be  computed  via  L 

independent  modulo  Pj^  operations.  A brief  overview  of  the 

Residue  Numbering  System  ( RNS ) is  provided  in  Appendix  A. 

The  set  of  complex  numbers  {a+b/T;  a,beZ}  form  a 

subring  generated  by  Z and  i=  /T  and  are  called  the  ring  of 

Gaussian  integers.  In  the  RNS,  the  corresponding  sets  of 

residues  modulo  pj^  are  defined  as  { a^^+bj^i  ; aj^,bj^eZp  }, 

k=l,2,...,L  and  each  forms  a field  (Z  [ i^ ],+,’, 0 , 1 ) under 

2 

multiplication  and  addition.  If  i^s-lmodpj^  has  no  solution 
then  ij^/Z  and  -1  is  called  a quadratic  non-residue 

K 

modulo  p.  . A numbering  system  for  which  i./Z  for  all  k 
K ^ Pk 

is  referred  to  as  the  complex-RNS  (CRNS).  In  the  CRNS  each 

aj^+bj^ij^,  k = l,2,...,L,  is  considered  to  comprise  real  and 

imaginary  parts  where  each  product  requires  four  real 

modulo  Pj^  multiplies  and  two  real  modulo  Pj^  adds. 

2 

Conversly,  if  ij^s-lmodpj^  has  a solution,  then  ij^eZ^j^  and  is 
a quadratic  root , -1  is  a quadratic 

residue . A numbering  for  which  ij^eZ^j^  for  all  k is 
referred  to  as  the  quadratic  RNS,  or  QRNS,  and  the  residue 
values  aj^+ij^bj^  are  treated  as  real. 
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If  ij^  is  a quadratic  root  then  there  exists  an 

2 2 

isomorphism  a:  [ i ] -»Zpj^  where  Zpj^  = Zpj^X  Z^j^  is  the  product 

set  of  Zpj^  with  itself.  Each  product  forms  a corresponding 

2 

set  of  two-tuples  ^pk~  ^ ^ ^ where  the 

rules  of  composition  satisfy 


( aj^  ,bj^ ) (|)(  Cj^,dj^  ) = ( ( aj^<J)Cj^ ) , ( bj^(|)dj^ ) ) ; k=l,2,...,L, 


for  <!>=(  + ,-,*),  each  taken  modulo  Pj^ . Here  addition  also 
requires  two  real  adds;  however , multiplication  requires 


only  two  real  multiplies  and  real  additions. 


The 


isomorphism  a;Zpj^->Zpj^  and  its  application  to  high-speed 


been 


covered 


concurrent  architectures  for  DSP  have 
extensively  in  the  literature  [55-57,60] 

If  Pj^  is  a Gaussian  prime  pj^=4K+l,  or  a composite  of 
Gaussian  primes,  the  isomorphism  is  the  following. 


• * 

a(  Uk'^'^k^k^  ^ ^k'^^k^k  ^ ®k"^k^k  ^ ^ '^k  ' ^k  ^ ' 

aj^=(  2"^  ( Zj^  + z^  ) )modpj^ ; bj^  = 2~^  j“^  ( Zj^-zj^ ) )modpj^  ; 

2 

jj^=-lmodpj^ ; k = l,2,...,L, 


(4-10) 
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where  is  a quadratic  residue.  Although  this  is  a 

very  recent  innovation,  a number  of  multi-moduli  papers 
(L>2)  have  appeared  in  print.  The  significant  advantage  to 
the  QRNS  is  the  computational  reduction  allowed  by  (4-10) 
and  the  term  QRNS  as  used  in  the  engineering  literature 
implies  application  of  the  isomorphism  in  (4-10). 

4.3.3  Single  Modulus  System 

Recently,  Taylor  [,58]  developed  the  QRNS  for  a single 
modulus  (L=l)  to  form  what  is  called  the  SM-QRNS . In  the 
SM  system  Taylor  was  able  to  overcome  the  principle 
performance  limitation  to  most  RNS  systems,  that  of  residue 
to  decimal  conversion.  The  moduli  used  to  define  the 

SM-QRNS  were  primes  of  the  form  p=2^+l  where  admissable 
values  of  n are  2,  4,  8,  16,  or  32.  In  this  form  the 
SM-QRNS  parameters  are  highly  composite  and  essentially 
radix-2  in  value  (i.e.  require  a simple  correction  stage). 

For  p=2'^+l  the  coefficients  in  the  isomorphism  (4-10) 
equate  to  2“^  = 2”~^+l;  and  ( 2 j ) “^=2"+l-2^/^  . Each 

of  these  SM-QRNS  scaling  modules  can  be  implemented  as 
binary  shift  and/or  operations  as  below: 

scaling  by  j 

jxn,odp=2"/2xj_g-x^j : 

scaling  by  2~^ 

2“^Xmodp=(  2^~^  + l )Xj^Q+X^j  ; X=2Xj^j+Xj^q 
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scaling  by  ( 2 j ) 


-1 


(2j)  ^modp=X^j-2^''/2  + 


HI  LO' 
(4-11) 


In  order  to  compete  with  traditional  complex 
arithmetic  units,  the  computational  units  that  make  up  the 
SM-QRNS  DWD  processor  must  be  designed  to  operate  with  high 
speed  and  low  hardware  complexity.  To  explore  this 
question,  a detailed  analysis  of  each  QRNS  single  modulus 
subsystem  architecture  was  performed  in  [58].  Based  on  an 
elementary  delay  of  a primitive  gate  (e.g.  NAND  TTL  gate), 
these  architectures  were  reported  to  integrate  into  an 
efficient  arithmetic  unit  with  a substantially  reduced 
hardware  latency  budget.  To  test  this  hypothesis  a gate 
array  SM-QRNS  unit,  based  on  GE's  IGC20000  technology  was 
undertaken  [61].  The  timing  studies  in  [61]  confirmed  the 
conjectures  advanced  in  [58].  Based  on  these  computational 
primatives,  a high  resolution  complex  intensive  multiply 
DWD  estimator  can  be  designed  to  achieve  high  throughputs. 
Furthermore,  the  symmetry  property  of  the  DWD  kernel  and 
the  corresponding  real-valued  nature  of  the  DWD  output  will 
be  shown  to  integrate  well  with  the  SM-QRNS. 
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4.2.4  Wigner  Processor  Implementation 

Figure  4-2  illustrates  the  SM-QRNS  DWD  processor.  The 
input  data  to  the  SM-QRNS  encoder  are  assumed  to  be  n-bit 
2's  complement  words.  Each  number  is  mapped  into  an 
(n+l)-bit  integer  modulo  p.  If  the  input  value  X<0,  then 
-X=TC(X)+2'^X  (TC  denotes  a 2's  (radix)  complement 
operation)  having  an  n-bit  format.  The  same  integer  in  the 
SM-QRNS  case  is  given  by  -Xs ( 2'^+1-X ) modp=TC ( X ) +1 . For 
example,  if  p=17  and  X=-l,  then  TC(1)-»1111  and 
TC(l)+l  = 16s-lmodl7.  If  N=-2,  then  TC(2)->1110  and 
TC( 2 )+l=llll=15s-2modl7 , etc.  Therefore,  the  code 
converter  can  be  implemented  by  modifying  a conventional 
2's  complement  encoder  with  a simple  increment  by  1 
ci rcui t . 

Data  windowing  is  performed  prior  to  two-tuple  QRNS 
conversion.  The  complex  windowed  data  points  x^(m)  and 
x^(-m)  are  mapped  into  their  respective  SM-QRNS  two-tuple 
equivalents  via  the  isomosphism  in  (4-10)  for  L=l. 
Computation  of  the  first  two-tuple  element,  z,  requires  two 
modular  operations:  j-magnitude  scaling  and  modulo  p 
addition.  The  first  operational  component,  the  j scalar, 
performs  scaling  on  b by  2^^^mod(  2'^+l ) which  is  implemented 
as  in  (4-11).  The  modulo  p adder  comprises  three  parts: 
an  n-bit  fast  carry-lookahead  adder,  a modulo  p mapping 


Figure  4-2 


SM-QRNS  DWD  processor  block 
diagram. 
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unit,  and  control  logic  (compare  network)  [58],  The  n-bit 
adder  is  a conventional  full  carry-lookahead  adder.  The 
compare  network  determines  whether  modulo  p mapping  is 
necessary.  If  S>p,  then  the  modulo  p mapping  unit  performs 
the  map  S<-S-p  using  the  procedure  detailed  in  reference 
[52]  : 

★ 

The  second  two-tuple  element,  z , requires  inclusion 
of  a negator  prior  to  modulo  p addition,  where  the  negator 
completes  the  following  operation:  x->-x=  ( p- | x | ) modp  . 

The  only  modular  components  necessary  to  generate  the 
DWD  kernel  are  modulo  p multipliers.  The  modulo  p 
multiplier  is  realized  via  an  unsigned  multiplier  and  a 
network  for  special  case  processing  [58].  For  inputs  x and 
y the  special  cases  are  the  following: 

if  both  X and  y equal  2’^,  then  set  their  product  to  1; 

if  X equals  2^  but  y does  not,  then  set  their  product 

to  -y, 

where  the  negator  is  used  to  obtain  -y.  If  the  single 
modulus  is  sufficiently  small  (p=2^+l;  n<8),  table-lookup 
methods  may  be  employed  where  a 2(n+l)-bit  address  is 
compressed  into  a 2n-bit  address  by  using  logic  gates  to 
discern  the  case  2^^  from  0.  In  [58]  Taylor  describes 
efficient  direct  lookup  methods  to  extend  the  allowable 
moduli  to  satisfy  n<16. 

Complex  conjugate  multiplication  is  completed  in  the 
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SM-QRNS  by  multiplying,  modulo  p,  the  first  two-tuple 
element  of  each  data  point  to  the  second  two-tuple  element 
of  the  other  data  point.  Each  operation  of  this  form 
produces  effectively  two  complex  kernel  output  values. 
That  is,  by  definition  the  DWD  is  conjugate  symmetric, 
which  in  the  SM-QRNS  translates  to 

(Zj.(n),z^  (n))  = (z^  ( -n  ) , z ^ ( -n  ) ) , therefore,  each  QRNS 
multiply  provides  two  kernel  output  values  where  the  kernel 
values  are  related  by  the  transposition  of  the  two-tuple 
elements.  The  appropriate  DWD  kernel  output  is  magnitude 
scaled  then  supplied  to  the  SM-QRNS  FFT  modules  for  each 
two-tuple  element.  After  the  first  complete  pass  of  data 
through  the  FFT  modules  DWD  kernel  generation  for  the  next 
consecutive  window  center  position  is  initiated. 

The  magnitude  scaling  module  (Figure  4-3)  comprises 

three  stages:  conversion  from  the  two-tuple  values; 

scaling;  then  conversion  back  into  the  two-tuple  values. 

Magnitude  scaling  is  performed  prior  to  each  butterfly 

stage.  Effectively  then,  magnitude  scaling  is  associated 

with  each  multiply  in  order  to  control  the  geometric 

dynamic  range  growth  of  the  RNS  products  (analogous  to  full 

precision  products  in  conventional  systems).  For  example, 

let  the  output  of  the  two  complex  butterfly  inputs  and 

Z-,,  where  Z.=X.  + iY.,  be  Z=A+iB  and  let  Z be  encoded  as  a 
2'  3 J 3 

signed  SM-QRNS  word  with  respect  to  a modulus  p=2'^+l 
(prime)  such  that  (p-l)/2  = 2'^ 


Then  set 


123 


n/2-1  n/2 

2 2 


n-1 

2 +1 


Figure  4-3. 


Magnitude  scaling  module. 
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max( |Xj I )=max( I Yj I )=V,  where  max { | A | ) =max ( | B | ) 

=2V^+Ng2V<2^  ^ for  Ng  the  number  of  additive  stages  in  the 
butterfly  module.  At  the  initial  pass,  multiplication  is 
by  one  and  magnitude  scaling  (associated  with  the  additive 
operations)  is  incorporated  into  the  second  level. 

The  SM-QRNS  FFT  modules  for  each  two-tuple  element  are 
independent  and  identical  by  definition  of  multiplication 
and  addition  in  (4-9).  Relative  to  a conventional  radix-4 
FFT  design,  for  example,  (Figure  4-4),  the  QRNS/FFT 
butterfly  reduces  from  a four  stage  process  (made  up  of  12 
real  multiplies,  followed  by  three  stages  of  add/subtracts) 
to  two  independent  three  stage  processes  (each  of  3 
real  multiplies  followed  by  four  add/subtracts  (+  one  -j 
scaling  operation)  then  four  add/subtracts )( Figure  4-5). 
The  SM-QRNS  design  relies  on  nearly  modulo  2*^  hardware  for 
all  derived  mathematical  operations  [54]. 

Recall,  the  DWD  is  a real-valued  function;  therefore, 
no  conversion  from  SM-QRNS  is  needed  at  the  output.  This 
fact  can  be  coupled  with  the  independent  and  identical 
relation  between  the  SM-QRNS  FFT  modules  for  each  two-tuple 
element  to  increase  the  throughput  in  the  DWD  processor. 

For  clarity,  consider  a power-of-2  FFT.  The  butterfly 
equations  for  x complex  express  as 


x^'^^d)  = x®(  i )+W*^x^(  j ) 
x^'^^(j)  = x®(  j)-W^x^(i) 


(4-12) 
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Four  point  FFT  module  for  mod2'^ 
based  processor. 


Figure  4-4 
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Four  point  FFT  module  for 
SM-QRNS  based  processor. 
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where  the  superscripts  denote  the  butterfly  stage  and 
the  appropriate  twiddle  factor.  In  an  SM-QRNS  two-tuple 
equivalent  (z,cz)  the  corresponding  butterfly  equations  are 
the  following 


.,5  + 1 , . . 


cz^^l(i) 


= cz® ( i )-cz^cz^( j ) 


™S+1 , . . 

"x  (3)  = 


cz^^l(j) 


= cz®( j )-cz^cz^( i ) 


(4-13) 


where  the  subscripts  denote  the  correspondence  of  the  QRNS 
value  to  its  integer  equivalent  (ie  x ->  (z^,cz^))  and  the 
superscripts  are  as  before.  Notice,  only  the  first 
two-tuple  element  z is  needed  at  the  input  stage  to  each 
butterfly  to  produce  the  output  z. 

For  real  output  the  first  element  of  the  two-tuple  in 
a SM  system  is  equivalent  to  the  integer  value;  therefore, 
only  the  first  two-tuple  element  needs  to  be  computed 
(naturally  z=cz)  and  because  the  DFT  operations  for  each 
two-tuple  element  are  independent,  only  those  associated 
with  the  first  two-tuple  element  are  needed  in  the  DWD 
processor  implementation  provided  DFT  processing  is 
performed  entirely  in  the  SM-QRNS.  That  is,  for  the  single 
modulus  system  magnitude  scaling  can  be  performed  via 
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table-lookups  applied  directly  to  the  first  two-tuple 
values,  eliminating  the  need  for  conversion  from  the 
SM-QRNS  two-tuple  elements  during  intermediate  processing. 
Here  the  SM-QRNS  butterfly  module  for  the  second  two-tuple 
element  is  employed  to  speed  up  processing  by  applying 
first  two-tuple  elements  and  supplying  their  appropriate 
SM-QRNS  twiddle  factors. 

Alternatively,  two  DWDs  can  be  simultaneously  computed 
using  one  SM-QRNS  DFT  by  composing  alternating  two-tuple 
elements  for  even  and  odd  values  of  t,  where  t denotes 
window  position,  such  that  each  two-tuple  applied  to  the 
DFT  is  of  the  form  t=0 , 2 , 4 , . . . , N-2  . The  first 
element  of  the  corresponding  two-tuples  at  the  output  of 
the  DFT  are  the  DWD  values  at  t and  the  second  two-tuple 
elements  are  the  DWD  values  at  t+1. 

For  longer  wordlengths  (n>16)  where  table-lookup 
scaling  operations  are  not  feasible,  conversion  from  the 
two-tuples  may  be  necessary.  Here  both  two-tuple  elements 
must  be  maintained  throughout  the  DFT  processor  until  the 
final  stage  of  processing  during  which  time  only  the  first 
two-tuple  element  is  needed.  At  the  final  pass  of  data, 
the  SM-QRNS  components  for  the  second  two-tuple  elements 
are  employed  to  speed  up  processing. 

As  the  DWD  avoids  the  traditional  STFT  trade-off 
problem  between  window  length  and  the  assumption  of 
stationarity , longer  observation  times  may  be  chosen  for 
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DWD  processing.  An  essential  part  of  any  RNS  based  system 
is  magnitude  scaling  and  is  important  to  the  choice  of  FFT. 
For  longer  transform  lengths  a large  number  of  nested 
multiplies,  associated  with  the  radix-2  FFT  and  radix-4 
structures,  translates  to  extensive  magnitude  scaling 
requirements,  which  should  be  avoided  in  an  RNS  system. 
Instead,  architectures  with  fewer  nested  multiplies  should 
be  considered.  Taylor  and  Huang  investigated  the 
comparison/trade-off  characteristics  associated  with 
various  standard  transform  alternatives  [53].  In  their 
work,  the  WFTA  proved  favorable  in  terms  of  reducing  nested 
scaling  operations;  however,  the  irregular  data  flow  and 
extension  requirements  are  a definite  disadvantage. 

A hybrid  system  which  uses  a large  radix  Cooley-Tukey 
decomposition  with  a small-WFTA  at  each  level  can  be  used 
to  reduce  nested  multiplies  while  avoiding  the  major 
di sadvantagees  associated  with  the  large-WFTA.  For 
example,  a large  N point  DWD  processor  using  a radix-L  CT 
decomposition  would  require  only  logj^N  nested  levels  of 
complex  multiplication  associated  with  the  twiddle  factors. 
Each  of  these  will  be  computed  in  the  SM-QRNS.  If  the 
radix-L  DFT  is  computed  via  an  L point  small-WFTA  one 
additional  multiply  stage,  consisting  of  L multiplies,  is 
added  to  each  level.  For  the  SM-QRNS,  the  j-scale 
operations  are  absorbed  into  the  multiply  stage  of  the 
small-WFTA  and  pre-weave  and  post-weave  operations  are 
exclusively  addition  and  subtraction. 
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4.3.5  Throughput  Comparison 


Assuming  all  worst 

case  conditions. 

for 

the  SM- 

QRNS 

(Table  4-1)  a direct 

comparison  was 

made 

between 

the 

SM-QRNS  based  DWD  processor  and  a 

conve 

ntional 

DWD 

processor  equivalent.  To  establish  a common  denominator, 
all  operations  were  treated  as  consecutive  and  both 
processors  were  assumed  to  have  input  data  in  a 2's 
complement  word  representation  and  fixed  integer  scaling. 
Tables  4-2  and  4-3  provide  respectively  the  throughput 
increase  for  the  SM-QRNS  DWD  processor  for  the  case  of  a 
radix-2  SM-QRNS  FFT  module  and  a radix-4  SM-QRNS  FFT 
module.  Percentage  throughput  is  defined  as 

SM-QRNS- (^^°/")mod2>^ 

n 

(DWD/A)mod2"  (4-14) 

where  A is  unit  of  delay.  The  wordlength  for  which  the 
SM-QRNS  becomes  faster  is  dependent  on  the  transform  length 
and  the  FFT  butterfly  size  and  is  n>5  or  less.  For  such 
small  word  sizes,  however,  table-lookups  may  be  used  to 
perform  magnitude  scaling  and  only  the  first  two-tuple 
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Table  4-1,  Analysis  of  SM-QRNS  DWD  processor 
throughput . 


Assumptions  in  Comparison 


n 

(1)  comparison  is  with  mod2  DWD  processor 

(2)  magnitude  scaling  is  required  for  each  multiply 

(3)  initial  conversion  and  conversion  from/to  SM-QRNS 
for  each  magnitude  scaling  operation  (i.e.  no 
table-lookups 

(4)  both  two— tuple  elements  must  be  maintained 
throughout  the  processor,  until  the  final  pass 
of  data  (this  follows  from  (3)) 

(5)  all  multiplies  are  negative  (negative  multiplies 
have  greatest  delay  in  SM-QRNS) 

(6)  the  symmetry  property  of  the  DWD  kernel  is 
applied  to  BOTH  processors  to  reduce  computations 

n 

(7)  in  the  mod2  processor,  the  delay  associated  with 
complex  conjugation  is  equivalent  to  that  of 
negation 

(8)  operations  are  treated  as  consecutive 
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Table  4-2.  The  percentage  increase  in  throughput 
for  the  SM-QRNS  given  a radix-2  FFT 
as  the  bases  for  the  DFT  stage  in 
both  processors. 


TRANSFORM 
LENGTH  = N 

WORDLENGTH 
= n 

N 

n = 16 

n = 32 

64 

54.6 

78.6 

256 

50.2 

74.4 

1024 

48.4 

71.9 

4096 

46.9 

70.3 

wordlength  at  which 

the  processing 

speeds 

become  equivalent: 

n < 5 
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Table  4-3.  The  percenatge  increase  in  throughput 
for  the  SM-QRNS  given  a radix-4  FFT 
as  the  bases  for  the  DFT  stage  in 
both  processors. 


TRANSFORM 

WORDLENGTH 

LENGTH  = N 

= n 

N 

n = 16 

n = 32 

64 

71.9 

97.1 

256 

62.9 

87.5 

1024 

57.8 

82 . 0 

4096 

54.5 

78.5 

wordlength  at  which 
become  equivalent: 

the  processing 

speeds 

N 

n 

64 

<3 

256 

<4 

1024 

<4 

4096 

<5 
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element  is  needed,  translating  to  a factor  of  two  in 
throughput  increase  for  the  SM-QRNS. 

To  achieve  very  high  data  rates  the  FFT  stage  of  the 
processor  may  be  implemented  using  the  butterfly  pipeline 
algorithms  developed  in  [62]  and  [63]  for  respective 
radix-2  and  radix-4  FFT  implementations.  A variable-length 
delay  commutator  design  approach  which  overcomes  most  v of 
the  complexity  problems  associated  the  pipeline  algorithms 
was  provided  in  [64]  and  used  to  develop  a delay  commutator 
circuit  for  a radix-4  FFT  pipeline  processor.  For  the 
radix-4  QRNS/FFT  processor,  SM-QRNS  two-tuple  elements  for 
each  set  of  4 DWD  kernel  values  are  pipelined  through  their 
respective  butterfly  modules,  a SM-QRNS  scaling  module,  and 
two  SM-QRNS  real  word  delay  modules. 

4.2.6  Weighted  Wigner  Kernel  Generation 

The  spDWD  complex  kernel  is  explicitly  complex 
computation  intensive.  The  high  speed  concurrent 
architecture  of  Figure  (4-2)  is  used  to  generate  the  spDWD 
kernel.  As  with  the  DWD,  only  one  half  the  spDWD  kernel  is 
computed.  The  values  for  the  other  half  of  the  spDWD 
kernel  are  determined  by  transposing  two-tuple  elements  at 
the  input  to  the  SM-QRNS  FFT. 

SM-QRNS  architectures  for  each  of  the  modular 
components  that  make  up  the  SM-QRNS  spDWD  kernel  processor 
have  been  provided  by  Taylor [58]  for  application  to  a SM 
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Complex  ALU.  Theoretically  produced  performance  measures 
were  conjectured  in  [58]  then  experimentally  confirmed  in 
[61].  Based  on  these  operation  times  in  a conventional 
radix-2  system  the  complex  multiplier  cells  would  complete 
a complex  multiplication  on  the  order  of  (28n-4)Tg,  for  Tg 
a unit  gate  delay,  and  the  multiply-add  cells  would  perform 
at  (14n+6)Tg.  In  the  SM-QRNS  the  time  to  complete  each 
complex  conjugate  multiply  will  be,  under  worst  case 
conditions,  (14n+34)Tg.  If  the  time  to  complete  conversion 
from  SM-QRNS  is  included  at  the  complex  conjugate 
multiplier  cells  the  total  processing  time  at  that  stage 
will  be  on  the  order  of,  again  under  worst  case  conditions, 
(14n+85)Tg.  Assuming  the  time  to  effect  complex 
conjugation  is  essentially  the  same  as  negation,  the  total 
processing  time  at  the  first  stage  of  the  radix-2  processor 
is  (28n+l)Tg. 

The  cells  do  not  operate  until  all  operands  have  been 
supplied;  therefore,  processing  speed  is  primarily 
dependent  on  the  speed  of  the  slowest  cell.  In  the  radix-2 
system  the  first  stage  of  the  smoothed  kernel  processor 
will  greatly  impede  system  throughput,  almost  by  a factor 
of  one-half,  whereas  inclusion  of  the  SM-QRNS  cells 
produces  like  processing  speeds  for  large  n,  minimizing  the 
amount  of  time  cells  are  standing  idle.  For  smaller 
wordlengths,  kernel  generation  is  performed  entirely  in  the 
SM-QRNS  where  magnitude  scaling  is  incorporated  directly 
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into  the  table-lookup  operations.  The  first  two-tuple 

element  for  the  entire  kernel  is  carried  through  the  kernel 

generation  process,  afterwhich  the  second  two-tuple  element 

★ 

is  obtained,  if  necessary,  via  the  relation  (n)=z^(-n). 

If  speed  is  a priority,  table-lookups  can  used  to 
perform  the  modulo  p multiplication.  Here,  magnitude 
scaling  can  be  incorporated  into  the  table-lookups.  If  the 
number  of  window  selections  is  reasonably  limited  (ie 
rectangular,  Hamming,  Blackman,  Gaussian),  table-lookups  in 
the  form  of  scaling  operations  can  be  used  in  place  of 
multiplications  in  the  multiply-add  cells,  effectively 
reducing  the  table  size  by  2 where  2 is  the  number  of 

available  window  selections. 

4 . 3 High-frequency,  High-resolution  Wigner  Processing 

As  with  the  continuous-time  case,  filtering  and 
modulation  are  important  DSP  operations.  Fortunately,  if 
the  respective  spectrums  of  both  x(n)  and  h(n)  are  assumed 
to  vanish  over  one-half  the  interval,  either  by 

oversampling  or  in  a case  as  analytic  signals,  the 
properties  of  filtering  and  modulation  are  analogous  to 
their  continuous-time  counterparts.  If  W (t,k),  W (t,k), 

y ^ 

and  Wj^(t,k)  are  the  respective  DWDs  of  y(n),  x(n),  and 


h ( n ) , then 
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y(n)=x(n)*h(n) 


N-1 

=»W  (n,k)=l/N[j;  W (n-m,k)W.  (m,k)  ]R.,(m)  (4-15) 

^ m=0  X n N 


and , 


y(n)=x(n)h(n) 


N-1 

=»W  (n,k)=l/N[  E W^(n,k-i)W,  (n,i)  ]R.,(i) 
y X n N 


where 


R 


N 


0<n>N-l 
otherwise . 


(4-16) 


The  filtering  and  modulation  properties  can  be  used  to 
produce  a band-select  Wigner  spectrum.  That  is,  like  the 
FFT  the  DWD  is  a baseband  analysis  tool.  Many  processing 
applications,  however,  involve  concentration  on  a selected 
band  of  frequencies  away  from  baseband.  To  achieve 
sufficient  spectral  resolution  in  such  cases  may  require 
more  sample  points  than  can  be  handled  efficiently  by  the 
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DWD  analyzer.  However,  because  an  FFT  is  imbedded  in  the 
DWD  stucture,  a band-selectable  (or  pseudo-zoom)  DWD  can  be 
designed  to  produce  fine  resolution  spectral  analysis  over 
a band  of  prespecified  frequencies  [65]. 

Figure  4-6  illustrates  the  pseudo-zoom  process.  Prior 
to  computation  of  the  Wigner  kernel,  the  data  is  bandpass 
filtered  over  the  frequency  band  of  interest;  this  so 
called  "observation  band"  is  heterodyned  to  baseband, 
lowpass  filtered  then  decimated  in  time.  The  DWD  of  the 
modulated  output  is  then  computed. 

The  bandpass  filter  passes  only  the  observation  band 
and  prevents  overlap  of  the  signal  into  the  observation 
band  during  the  heterodyning  process.  The  DWD  is  a 
real-valued  function  which  does  not  preserve  phase 
information  of  the  signal.  That  is,  recovery  of  the  signal 
from  the  DWD  can  be  obtained  only  up  to  a constant.  For 
this  reason  if  the  signal  type  to  which  the  process  is 
being  applied  is  unlikely  to  undergo  a 180  degree  change  in 
phase  over  the  observation  band  relaxed  phase  conditions 
may  be  warrented  to  improve  the  passband  flatness,  and 
rolloff  and  attenuation  characteristics.  Otherwise,  linear 
phase  filters  are  needed. 

The  term  "pseudo"  is  used  because  by  filtering  out  all 
but  the  frequency  band  of  interest  some  cross-terms 
resulting  from  the  inner  product  that  would  normally  fall 
within  the  observation  band  have  been  eliminated.  Recall 
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x(  n) 


Figure  4-6.  Block  diagram  of  band-select 
DWD  process. 


>W( n, k ) 
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(Equation  (3-22)),  when  the  Wigner  )<ernel  is  formed  all  of 
the  sinusoidal  components  are  mapped  into  the  form 
l/2(cos2w+l)  where  co  is  the  original  frequency  of  the 
sinusoid,  or  into  the  form  cos  ( 00^+002 ) t , where  00^  and  (*>2  are 
the  frequencies  of  two  component  sinusoids  with  two 
different  frequencies.  If  is  within  the  observation 
band,  however  0^2  is  not,  the  cross  term  contributions  at 
(Oj^+o)2  do  not  appear  in  the  DWD  of  the  band-select  signal. 

For  real  signals,  frequencies  are  mapped  to  dc  during 
computation  of  the  Wigner  Icernel,  effectively  corrupting 
the  signal  around  dc,  and  dc  information  is  difficult  to 
obtain  from  the  DWD  (Section  3.3).  Therefore,  insertion  of 
a "guard  band"  is  needed  during  the  heterodyning  process; 
that  is,  the  observation  band  is  modulated  to  "near" 
baseband,  where  the  guard  band  is  system  and  application 
dependent  [65].  For  example,  consider  the  case  of 
rectangular  windowing  with  acceptable  spectral  degradation 
-26db:  for  the  rectangular  window  an  expected  rollof  rate 
is  13db  per  sidelobe  and  the  guard  band  must  be  greater 
than  the  width  of  the  second  sidelobe.  For  analytic 
signals,  the  guard  band  is  not  needed  (Section  3.3). 

The  lowpass  filter  passes  only  frequencies  within  the 
analysis  band  plus  the  guard  band  and  filters  out  any 
unwanted  frequencies  outside  the  band  of  interest  that 
result  from  folding  during  the  heterodyning  process.  As 
with  the  bandpass  filter,  the  filter  characteristics  are 
system  and  application  dependent. 
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Finally,  the  deciraator  is  used  to  reduce  the  effective 
sampling  rate  of  the  input  to  the  DWD  processor  to  achieve 
increased  resolution.  The  decimator  type  is  also  system 
and  application  dependent  and  is  primarily  a trade-off 
between  resolution  efficiency  and  system  complexity. 

After  the  band-select  process  is  complete,  the 
preprocessed  output  is  input  to  the  kernel  generator  of  the 
DWD. 

To  illustrate  the  "pseudo"-zoom  process,  a music 
signal  sampled  at  20kHz  (Figure  4-7)  will  be  observed  over 
a IkHz  band  between  5500  and  6500Hz.  The  filter 
characteristics  used  in  the  band-select  DWD  example  are 
given  in  Figures  4-8  and  4-9.  Decimation  is  effected  as  a 
simple  switch  and  performs  a decimation  ratio  of  4:1. 
Figure  4-10  is  a 64  point  DWD  of  the  original  music  data 
and  Figure  4-11  is  the  corresponding  band-select  64  point 
DWD  output  for  the  observation  band,  5500  to  6500Hz. 
Spectral  resolution  over  the  observation  band  is  not 
obtained  for  a 64  point  DWD,  however  the  "pseudo"-zoom 
process  resolves  the  individual  tones.  For  comparative 
purposes  a 256  point  DWD  of  the  music  data  and  a 256  point 
DWD,  filtered  over  the  observation  band,  of  the  music  data 
are  given,  respectfully,  in  Figures  4-12  and  4-13. 
Observe,  in  the  256  point  DWD,  the  presence  of  tones  that 
are  not  present  in  the  band-select  64  point  DWD.  These 


142 


VALUES 


Figure  4-7.  Music  signal  sampled  at  20kHz. 
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BNO  PASS  ELLIPTIC  FILTER 

PASS  BAND  RIPPLE  0.200  08 

PASS  8AND  EDGES  5S00.001  S500.000  KZ 

STOP  BAND  EDGES  3092.771  £833.496  HZ 

SAMPLING  frequency  20000.000  KZ 

ATTENUATION  -40.000  08 


OENOMimTOR 

1 1 .0OOOOOE4-OO 

2 2.323302E-F00 

3 3.470944E+00 

* 6.83S298E+00 

3 3.228624E+00 

6 6.093274E+00 

7 4.323278E+00 

8 1.629330E+00 

9 6.228603E-01 


NUMERATOR 
1 .233341E-02 
2.466397E-02 
4.917343E-02 
3.984104E-02 
7.337381E-02 
3.984103E-02 
4.917343E-02 
2.466397E-02 
1 .233341E-02 


(OBI 


Figure  4-8.  Bandpass  filter  for  "psuedo"-zoom 
DWD  example. 


SF84PLING  FREQUENCY 
LOW  PASS  CHEBYCHEV 
STOP  BAND  ATTENUAT 
BM40  EDGES 

DENOMINATOR 

1 l.OOOOOOE+00 

2 -S.313048E+00 

3 1.739896E+01 

4 -2.824037E+01 

3 2.831816E+01 

6 -I .834670E+01 

7 7.382170E+00 

8 -I .7807800+00 

9 1.839127E-01 


20000.000  HZ 
FILTER 

ON  40.000  DB 

1100.000 

numerator 

9.093963E-03 
-3.067343E-02 
1 .337898E-01 
-2.22S362E-01 
2.607683E-01 
-2.223362E-01 
1.337897E-01 
-5.067343E-02 
9.095962E-03 


(Utl) 


Lowpass  filter 
DWD  example. 


Figure  4-9 


for 


psuedo"-zoom 
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Figure  4-10. 


DWD  of  music  data 


for  N=64. 


fWGNITUDE 


Figure  4-11. 


Band-select  DWD  of  music  data 
for  N=64. 
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Figure  4-12. 


DWD  of  music  data 


for  N=256. 


Figure  4-13. 


DWD  of  music  data,  filtered  over 
the  observation  band,  for  N=256. 
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MAGNITLCe 


Figure  4-14. 


Band-select  DWD  of  music  data 
for  N=256. 
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result  from  the  cross-terms  contributions  between  tones 
outside  the  observation  band.  Figure  4-14  is  a 256  point 
band-select  DWD  where,  as  expected,  the  cross-term 
contributions  have  been  elliminated. 

To  increase  resolution  and  reduce  computations  the 
simplified  band-select  DWD  of  Figure  4-6  may  be  modified  as 
follows : 

partial  decimation  =>  improved  attenuation  of  BPF; 

lowpass  filter  (complex  coefficients  =>  SM-QRNS 

candidate ) ; 

decimation ; 

complex  demodulation  to  "near"  baseband; 

kernel  generation; 

FFT. 

Working  with  linear  phase  filters  makes  achieving  a 
flat  passband  with  sharp  attenuation  a difficult  task; 
especially  if  the  observation  is  very  narrow  relative  to 
the  total  band.  To  prevent  aliasing  in  the  the  DWD,  real 
data  sequences  require  that  the  signal  be  sampled  at  double 
the  Nyquist  rate  and  analytic  signals  require  sampling  to 
be  at  the  Nyquist  rate.  For  bandpass  signals  then  aliasing 
can  be  prevented  with  a sampling  frequency  of  four  or  eight 
times  the  bandwidth  over  which  the  spectrum  is  nonzero. 
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(depending  on  whether  the  signal  is  real  or  analytic).  As 
such,  for  narrow  observation  bands  away  from  baseband, 
partial  decimation  can  be  performed  immediatly  prior  to 
bandpass  filtering.  By  doing  this  high  attenuation  of  the 
bandpass  filter  can  be  effected  for  a flat  passband.  For 
example,  consider  an  observation  band  of  400Hz  over  a total 
band  of  4KHz  (as  applicable  to  speech).  For  real  data  the 
sampling  rate  would  exceed  16KHz.  If  a 4:1  decimation 
where  performed  after  bandpass  filtering,  effecting  a 4KHz 
sampling  rate, the  bandpass  criteria  to  prevent  aliasing 
would  still  be  satisfied.  After  lowpass  filtering  then  a 
2:1  decimation  may  be  performed  producing  an  overall  8:1 
decimation.  Due  to  the  increased  attenuation  of  the 
bandpass  filter,  the  required  DWD  "guard-band"  is  reduced; 
this  translates  to  higher  resolution  in  the  "psuedo"-zoom 
DWD  output. 

If  complex  demodulation  is  used  on  the  signal,  the 
mixing  process  can  be  passed  through  an  FIR  lowpass  filter 
as  illustrated  in  Figure  4-15.  Here  the  coefficients  of 
the  LPF  are  complex  however  decimation  of  the  sequence  can 
be  performed  prior  to  the  mixing  process,  hence  reducing 
the  number  of  multiplies  in  the  mixing  process  by  the 
decimation  ratio.  Either  the  lowpass  filter  or  the  entire 


complex  demodulation  process  may  be  implemented  in  the 
QRNS.  In  either  case  conversion  to  QRNS  would  be  performed 
after  bandpass  filtering  and  partial  decimation. 
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Figure  4-15.  Multiply  reduction  a)Complex 

demodulation  followed  by  lowpass 
filtering;  b)Lowpass  filtering 
with  complex  coefficients 
followed  by  complex  demodulation. 
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4 . 4 Acousto-optic  Wigner  Processor 

Many  problems  requiring  enhanced  time-frequency 
analysis  require  computational  bandwidths  well  beyond  those 
which  can  be  obtained  using  conventional  digital 
implementations.  For  example,  electronic  countermeasure 
(ECM)  processors  and  a high  bandwidth  doppler  radar,  using 
a chirp  envelope,  fall  into  this  category.  To  intercept 
these  potential  problems,  an  alternative  technology  must  be 
considered.  In  this  section  an  acousto-optic  Wigner 
processor  configuration  is  presented  as  a potential 
real-time  signal  processor. 

Several  papers  have  suggested  producing  a Wigner 
spectrum  by  optical  means.  Bartelt  et.al.  [32],  in  1980 
presented  preliminary  experimental  results  on  an  optical 
processor  whose  output  is  related  to  an  approximation  of 
the  Wigner  spectrum;  it  was  called  a "local  spectrum." 
Their  embodiment  did  not  provide  a means  of  producing  the 
Wigner  kernel.  Kumar  et  al . performed  a paper  study  of  an 
optical  Wigner  processor  [66].  They  assumed  the  existence 
of  an  optical  Wigner  processor  and  conjectured  its 
performance  based  on  digital  simulation.  The  simulator  did 
not  emulate  the  optical  subsystem  but  simply  signal  in 
noise  tests  on  the  WD , This  section  will  address  more 
definative  optical  processor  design  issues.  More 
specifically,  in  this  section  an  acousto-optic  specialized 
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Wigner  distribution  processor  for  operation  at  very  high 
data  rates  will  be  developed. 

Acousto-optic  processing  techniques  have  been 
successfully  applied  to  spectrum  analysis  [67-69].  A laser 
beam  provides  a coherent  source  which  is  expanded  and 
recollimated.  The  expanded  beam  passes  through  a surface 
acoustic  wave  device  where  the  near  surface  interaction 
with  a diffraction  function  produces  the  transmittance  beam 
eminating  from  the  SAW  device. 

To  produce  the  diffraction  pattern,  a signal  s(t)  is 
applied  to  a transducer  at  one  end  of  the  SAW  device.  The 
transducer  generates  a surface  acoustic  wave  with 
propagation  delay  t=x/v,  where  x is  in  the  direction  of 
propagation  and  v is  the  velocity  of  the  acoustic  wave 
[68-71].  The  surface  acoustic  wave  travelling  along  the 
cell  effects  a moving  diffraction  pattern  of  the  form 
s ( t-x/v ) . 

Given  a laser  source,  the  incident  field  assumes  a 
Gaussian  plane  wave  of  the  form 


g(x)=a(x)e  ^“l^; 


( 4-17 


a(x)=/57i(l/aQ)e  ^ ^^0 


where  oo^  is  the  photon  angular  frequency  c/X  and  aQ  is  the 

The  acousto-optic  interaction  produces  a 


beam  waist. 
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refracted  wave  which  when  operating  in  the  Bragg  regime  can 
be  represented  by  a principle  diffracted  beam  in  quadrature 
with  an  undiffracted  main  beam  [68-71].  That  is,  the 
transmittance  field  can  be  approximated  as  a phase 
modulation  on  the  incident  optical  beam  to  produce  a 
complex  field  amplitude  of  the  following  form 


e ( t , X ) =a ( X 


jas( t-x/v) 


(4-18) 


where  a is  an  attenuation  term.  If 
small  in  magnitude,  the  diffracted 
by  a first  order  Taylor  series 
exponent  producing  the  transmittance 


the  acousti 
beam  can  be 
approximati 
field  ampli 


c wave  is 
represented 
on  on  the 
tude 


e ( t , X ) =a ( X ) [ 1 + j as ( t-x/v ) ] 


(4-19) 


which  comprises  the  diffracted  and  undiffracted 
components . 

To  achieve  optimum  diffraction  efficiency 
magnitude  of  the  photon  momentum  must  be  conserved, 
mismatch  of  the  incident  beam  and  the  transmittance 
minimized  when  the  optical  wave  is  incident  on  the 
acoustic  wave  at  the  Bragg  angle,  0g^=sin  ^(X/2nA), 
the  acoustic  wavelength  and  n the  refractive  index  o 


field 

, the 
Photon 
wave  is 
surface 
for  A 
f the 


153 


Figure  4-16. 


Bragg  cell  diffraction  of 
counterpropagating  waves. 
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cell  material  [68—71]  (Figure  4—16).  Given  the  acoustic 
wave  propagation  to  be  parallel  to  the  cell  edge,  for  an 
angle  of  incidence  on  the  cell  at  the  Bragg  angle 
0B=sin  ^(X/2A),  then  according  to  Snell's  law,  the  optical 
wave  is  incident  on  the  acoustic  wave  in  the  material  at 
the  Bragg  angle  in  the  device  material,  given  by  9 . A 
portion  of  the  light  is  diffracted  from  the  incident  beam 
at  an  angle  29^^.  Again  via  Snell's  law,  the  angle  between 
the  transmittance  beams  is  29„. 

The  Wigner  inner  product  is  produced  by  imposing  a 
double  diffraction  on  the  optical  beam  with  two 
counterpropagating  surface  acoustic  waves  referenced  from 
the  cell  center.  Figure  4-17  illustrates  the  acousto-optic 
Wigner  processor  configuration.  The  light  source  is 
incident  on  Cl  producing  the  diffracted  wave  complex 
amplitude 


e(t,x)=a'(x)s( t-x/v) 


(4-20) 


diffracted  at  an  angle  with  respect  to  the 
main  beam  of  29^.  The  weighting  term  a'(x 
truncated  form  of  a(x)  dependent  primarily 
aperture,  beam  expansion  characteristic 


undiffracted 
is  a scaled, 
on  the  cell 
and  acoustic 


attenuation  [68-71]. 
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EXPANDER  CELL 


Figure  4-17.  Optical  Wigner  processor 
configuration . 
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Spatial  filtering  to  eliminate  the  undiffracted  main 
beam  is  effected  by  L3,  Si  and  L4 . Lens  L3  images  the 
transform  of  e(t,x)  onto  Si.  As  e(t,x)  is  diffracted  at 
the  angle  20^  from  the  undiffracted  beam,  where  0^  is  given 
above,  they  are  spatially  separated  in  the  transform  plane 
by  an  amount  fX/A,  A slit  at  SI  passes  only  the  diffracted 
beam.  Lens  L4  images  the  undiffracted  beam  onto  cell  C2 

where  a double  diffraction  is  imposed  by  the 

counte rpropagat ing  acoustic  wave  at  the  Bragg  angle.  The 
transmittance  field  amplitude  comprises  a principle  beam  of 
the  form 

ee( t ,x)=a ' ( x ) s ( t-x/v ) s ( t+x/v ) , ( 4-21 ) 

diffracted  from  the  main  beam  by  20g.  An  optical  stop  at 
Dl  is  used  to  eliminate  the  undiffracted  term.  The  complex 
amplitude  at  the  back  focal  plane  of  the  transform  lens, 
L5,  expresses  as  [refer  to  Appendix  B for  details] 

00 

W(  t , Xj  ) =Kja '(  X ) s ( t-x/v ) s ( t+x/v ) e ^ ^ dx . (4-22) 


The  spatial  resolution  in  the  transform  plane  is  seen  to  be 
primarily  dependent  on  the  weighting  term,  a' (x) . If  the 
distance  transversed  by  the  collimated  beam  is  sufficiently 
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small  relative  to  the  length  of  the  cell  and  the 
collimating  lens  aperture,  the  smoothing  factor  a'(x) 
approximates  a constant  over  the  length  of  the  acoustic 
cell  [72]. 

For  analytic  signals  complex  conjugation  of  the  signal 
can  be  effected  via  quadrature  modulation  then  SSB 
filtering  [68].  For  s^{t,x)  and  Sj(t,x)  the  respective 
real  and  imaginary  signal  components,  the  quadrature 
modulated  signal  expresses  as 


s(t,x)  = s„  ( t , X ) cosco  x-s^  ( t ,x  ) sinco  X . 


(4-23) 


The  transform  of  s(t,x),  taken  with  respect  to  x,  then 
expresses  as 

[ Sj^(  ) +Sj^(  ) } - j { S j ) -S  j ( 03+0)^  ) } . ( 4-24  ) 

Equation  (4-24)  can  be  seen  as  simply  the  transform  of  the 

demodulated  signal  around  and  the  transform  of  the 

complex  conjugate  around  Placing  a screen  at  the 

focal  point  of  the  transform  lens  that  passes  only  the  term 

at  +w  then  inverse  transforming  the  wave  produces 
c 

s* ( t , X ) e”^“c^ . The  phase  term  isomorphically  translates 
the  output  spectra  an  amount  equivalent  to  in  the  Wigner 
plane  which  can  be  accomodated  for  at  the  detector. 


158 


Figure  4-18  indicates  a hybrid  optical  guided-wave 
Wigner  processor  configuration.  A planar  waveguide  is 
employed  where  surface  acoustic  waves  are  excited  on  a 
substrate  by  transducers  placed  at  opposite  ends  of  the 
device.  Acoustic  absorbers  are  placed  between  the 
transducer  and  the  device  edge  to  eliminate  unwanted 
counte rpropagating  waves  that  result  from  reflections 
[68-71].  Geodesic  lenses  and  Bragg  modulators  have 
successfully  been  fabricated  in  Ti-indif fused  LiNbO^ 
substrates  [69].  GaAs  substrates  have  received 
considerable  attention  toward  the  development  of  totally 
integrated  optical  processors;  however,  LiNbO^  substrates 
maintain  fewer  loss  characteristics  than  GaAs  substrates 
[68,69].  The  laser  diode  is  butt-coupled  to  the  waveguide 
edge  and  the  input  is  coupled  into  the  guide  which  forms 
the  transmission  medium.  The  major  component  problem  is 
one  of  two-dimensional  detection.  To  produce  an  output 
pattern  as  a function  of  t and  x^,  a parallel  readout 
linear  detector  array  [68]  can  be  coupled  to  the  waveguide 
edge  via  fiber  optics. 

In  the  above  optical  guided-wave  Wigner  processor 
configuration,  signals  of  small  BT  ( B=bandwidth , T=time) 
values  were  assumed,  as  the  proposed  processor  employs 
space-integrating  techniques  and  is  therefore  limited  by 
the  delay  line  length,  lens  aperture,  etc.  For  signals  of 
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Figure  4-18. 


Hybrid  optical  guided-wave 
Wigner  processor. 
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larger  BT  values,  a more  sophisticated  processor  employing 
time-integrating  techniques  is  necessary;  otherwise,  the 
pseudo(or  windowed)-WD  must  be  observed,  where  the  output 
is  a frequency  smoothed  version  of  the  WD . With  the  pWD, 
all  time  domain  properties  are  maintained.  Still,  the  pWD 
proves  an  attractive  alternative  to  STFT  techniques 
( Chapter  2 ) . 


CHAPTER  5 

SUMMARY  AND  CONCLUSIONS 


A time-dependent  spectral  estimator,  called  the  Wigner 
distribution  (WD),  has  been  studied.  Its  relation  to  the 
corresponding  one-dimensional  time  and  frequency  domains  of 
a signal  was  addressed.  The  effects  of  filtering  and 
modulation  on  the  WD  were  investigated  and  shown  to 
overcome  some  of  the  problems  associated  with  the  WD.  A 
Gaussian  weighting,  for  both  the  filter  and  window,  were 
shown  to  achieve  minimum  variance  in  the  weighted-WD  plane. 
The  relative  constraint  on  the  variance  between  window  and 
filter  required  to  maintain  positivity  was  determined  to  be 
for  a normal  distribution. 

The  double  Fourier  transform  relationship  between  the 
WD  and  the  AF  was  addressed.  Accordingly  the  WD  and  AF 
both  can  be  used  for  signal  interpretation  in  much  the  same 
way  the  one-dimensional  time  and  frequency  domains  of  a 
signal  are  both  used.  For  example,  the  duality  between  the 
AF  and  the  WD  can  be  applied  to  analysis  of  incomplete  data 
from  a nonstationary  signal  in  a manner  consistent  with 
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that  used  between  the  one-dimensional  time  and  frequency 
domains • 

The  DWD  is  an  FFT  based  estimator  and  is  therefore 
recognized  as  a viable  competator  to  the  STP.  A detailed 
comparison  between  the  DWD  and  the  STP  has  been  provided 
both  analytically  and  experimentally.  For  comparative 
purposes,  broadband  communication  signals  were  selected 
because  of  their  conceptual  accessibility  and  because  their 
dispersion  characterstics  are  highly  time  dependent. 

The  spDWD  was  interpreted  in  terms  of  filtering  and 
modulation  then  used  to  improve  DWD  estimation  of  signal 
types  for  which  the  DWD  produces  poor  results.  Again,  the 
signal  examples  were  broadband  communication  signals. 
Application  of  the  proposed  window  and  filter  shape  were 
experimentally  tested,  then  applied  to  a practical  example. 
Variations  of  the  filter  variance  for  a window  was  shown  to 
influence  the  estimation. 

In  practical  applications  the  signals  used  in  the 
examples  would  operate  at  higher  frequencies  than  used  in 
the  simulations.  The  reason  for  selecting  lower 
frequencies  was  primarily  for  illustrative  clarity.  That 
is,  if  the  signals  were  carried  at  the  higher  frequencies, 
for  64  point  and  256  point  transforms  a large  number  of 
consecutive  time  slices  would  be  required  to  display  one 
complete  period  of  the  signal.  This  would  serve  only  to 
obscure  the  display  without  contributing  to  the  analysis. 
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The  programs  used  to  effect  the  DWD  and  spDWD  were  written 
to  operate  interactively  with  the  ILS  subroutine  SDI . 

The  DWD  was  shown  to  require  sampling  be  at  double  the 
Nyquist  rate.  Therefore,  to  extract  equivalent  t-w 
information,  the  displacement  between  window  center 
positions  in  the  DWD  must  be  twice  that  of  the  STP.  The 
same  number  of  transforms  are  required  for  each,  however 
the  symmetry  of  the  DWD  kernel  was  shown  to  allow  a 
reduction  in  each  computation  as  will  be  discussed 
presently. 

Several  means  of  implementing  the  DWD  and  spDWD  for 
various  signal  types  were  presented.  For  real  input  data 
the  symmetry  of  both  the  DWD  kernel  and  the  corresponding 
DWD  output  kernel  was  shown  to  allow  computation  of  two  N 
point  DWD  time  slices  via  one  N/2  point  FFT.  This  was 
accomplished  by  rearranging  the  kernel  input  to  the  FFT  and 
applying  the  folding  properties  to  the  output.  Using  this 
result  and  modifying  the  approach  slightly,  computation  of 
each  DWD  slice  via  an  N/4  FFT  seems  to  naturally  follow. 
The  spDWD  kernel  is  also  real  symmetric  and  the  same 
results  apply  to  the  spDWD. 

For  complex  data  processing,  modular  arithmetic  was 
applied  to  the  processor  design  to  reduce  computations.  A 
numbering  system  called  the  quadratic  residue  numbering 
system  (QRNS)  was  chosen  because  it  allows  an  isomorphic 
mapping  from  the  complex  integer  representation  to  the 
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representation  of  two-tuples  where  multiplication  is 
completed  via  only  two  real  multiplies  and  real  one 
addition.  This  is  in  comparison  to  a standard  complex 
multiply  requiring  four  real  multiplies  and  two  real 
additions.  The  QRNS  for  the  case  of  a single  modulus  was 
used  in  the  design  where  the  chosen  modulus  was  of  the  form 
2^+1.  With  this  choice  of  modulus  the  isomorphism  to  the 
two-tuple  representation  was  a simple  set  of  shift/add 
operations.  Application  of  the  single  modulus  QRNS  was 
chosen  such  that  the  popular  radix-2  and  radix-4  FFT 
algorithms  could  be  used  in  the  processor  implementation. 
That  is,  in  a multi-moduli  system,  algorithms  requiring  a 
large  number  of  nested  multiplies  require  too  much  overhead 
associated  with  magnitude  scaling.  Naturally,  even  for  the 
SM  case  some  overhead  is  required  but  is  minimal.  A 
throughput  comparison  between  the  SM-QRNS  and  a 

conventional  system  showed  the  SM-QRNS,  for  a like  amount 
of  hardware,  to  have  higher  throughput  even  under  worst 
case  SM-QRNS  conditions. 

A band-select  DWD  for  high  frequency  analysis  was 
presented.  The  corresponding  elimination  of  cross-term 
components  was  discussed.  That  is,  cross-terms  that  result 
from  frequencies  outside  the  observation  band  did  not 
appear  in  the  band-select  implementation.  This  result 


should  be  accounted  for  when  interpreting  the  band-select 
output.  To  interpret  the  band-select  output,  comparison 
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with  the  spDWD  (for  the  spDWD  filter  bandwidth  less  than 
the  observation  band)  may  be  more  appropriate.  Recall, 
however,  the  spDWD  effects  a moving  filter  centered  at  each 
point  (t,w)  whereas  in  the  band-select  implementation  the 
filter  remains  centered  about  a fixed  frequency. 

Finally,  an  acousto-optic  processor  to  effect  the  WD 
was  presented.  The  processor  implementation  employed 
space-integrating  techniques  which  contrains  the  input 
signal  to  those  with  small  BT  values;  otherwise  the  pWD  is 
observed.  A processor  using  time-integrating  techniques  is 
needed  for  signals  of  larger  BT  values;  however,  a 
time-integrating  processor  implementation  is  not  as 
staightf orward , particularly  for  the  case  of  the  WD.  A 
possible  approach  may  be  to  effect  the  AF  (which  is  much 
simpler  for  the  case  of  time  integration),  then  use  the 
double  Fourier  transform  relationship  between  the  WD  and 
the  AF  to  produce  the  WD.  A processor  of  this  type  would 
most  likely  employ  both  time  and  space-integrating 
techniques . 


time 


APPENDIX  A 

RESIDUE  NUMBERING  SYSTEM  OVERVIEW 

The  definitions  and  relations  presented  here  are 
fundamental  to  abstract  algebra  and  can  be  found  in  a first 
year  graduate  text  in  the  area.  A particularly  thorough 

treatment  can  be  found  in  [73]  and  the  notation  used 

throughout  this  appendix  has  been  made  consistent  with  that 
used  in  [73]  for  textual  cross-referencing. 

The  integer  numbers  form  a monoid  over  both 

multiplication  and  addition  where  a monoid  is  defined 
as  a triple  (S,<j>,u)  in  which  S denotes  a nonvacuous  set,  <|> 
an  associative  binary  composition  in  S,  and  u a unique  unit 
element  for  which  a<}>u=a=u<t)a  for  all  aeS.  With  respect  to 
addition,  4)=+  and  u=0  while  for  multiplication,  <f>=*  and 

u=l;  for  both,  S=Z.  Any  subset  of  the  set  Z that  contains 

the  unit  element  and  is  closed  under  the  composition  (|) 
(i.e.  for  all  a,beS,  a<|)b£S)  is  referred  to  as  a submonoid 
( or  monoid ) . 

The  residue  numbering  system  (RNS)  representation  for 
integers  is  based  on  particular  equivalence  relations 
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termed  congruences  where  an  equivalence  relation  is  any 
relation,  E,  on  a set,  S,  which  satisfies  the  following 
properties : 

reflexive  (aEa;  aeS) 

symmetry  ( aEb  =>  bEa;  a,bcS) 

transitivity  (aEb  and  bEc  =>  aEc;  a,b,ceS). 

Now,  given  a monoid  (S,<J>,u),  a congruence , s,  in  S is 
an  equivalence  relation  in  S such  that  for  any  a,b,c,deS, 
if  a=b  and  csd  then  a<j)csb((>d.  Two  integers  are  defined  as 
congruent  modulo  a third,  denoted  a2b(modk),  if  a-b=nk 
where  n is  any  integer. 

The  definition  for  congruence  can  be  applied  to  the 
monoids  (Z,+,0)  and  (Z,*,l)  where  = is  ^(modk).  That  is, 
for  any  integer,  a,  there  exists  integers  b,k,aj^  such  that 
a=bk+aj^  for  k>0  and  0<aj^<k,  which  means  a=aj^(modk)  where  aj^ 
defines  a residue  class  (modulo  k)  given  by  ij^={aj^+nk; 
neZ}.  The  set  of  residue  classes  (modulo  k)  form  a 
subgroup  (or  group)  under  <|)  (where  <j)  can  be  used  to  denote 
multiplication  or  addition). 

A group  is  a monoid  in  which  all  of  its  elements 
are  invertible;  that  is,  for  every  aeS  there  exists  an 
inverse,  call  it  a such  that  aa  ^=u=a  ^a  where  a~^eS  and 
is  unique.  A subgroup  is  a submonoid  which  acts  as  a 


group . 
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The  group,  called  the  quotient  (or  factor)  group  of  Z 
relative  to  the  s(modlc)  and  denoted  Z/kZ,  is  of  order  k 
and,  hence,  is  isomorphic  to  the  group  under  <j)  of 
Zj^={0,l,  . . . ,k— 1}  . 

Two  monoids  (groups),  ( , <|)^  ,Uj^ ) and  ( $2  / <(>2  / ^2  ) » 

isomorphic  if  there  exists  a bijective  (i.e.  1 to  1, 

onto)  map  a:Sj^-»S2  such  that  for  a,beSj^ 

a(  a<(>j^b ) =a(  a ) <j>2a( b ) £$2  . The  map  a is  called  an  isomorphism 

of  onto  S2. 

By  definition  then,  a(u^)=U2  for  a an  isomorphism. 
The  isomorphism  a:Z/kZ-»Zj^  is  given  by  a(ij^)=aj^  where  each 
aj^eZj^  is  termed  a residue  and  equates  to  aj^  = amodk  for  each 
aeZ.  The  set  of  residues  is  loosely  termed  the  additive 
(multiplicative)  group  of  residue  classes  of  Zmodulok  and 
denoted  Zj^.  Although  informal,  this  notation  is  common  and 
will  be  used  as  given,  unless  stated  otherwise,  in  the 
remaining  discussion. 

Given  a finite  set  of  integers  { 0 , 1 , 2 , . . . , M-1 } = Z^^, 

each  integer  aeZ^^  has  a unique  RNS  representation  given  by: 

a“^(a^,a2/...,a^)j  aj^—amodpj^eSj^  (A  1) 


where 


^ ®k 

M=n  p,  ; gcd(p.  ,p  . )=1 . 

k=l  ^ ^ 


(A-2) 
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Signed  numbers,  ae [ -M/2 ,M/2 ) , are  represented  as  above  if 
a>0;  if  a<0,  then  aj^=  ( M- | a | ) modpj^ . Each  residue  set  Z^, 
k = l,2,...,L  forms  a subset  of  where  each  triple  (Zj^,  + ,0) 
defines  a subgroup  and  (Zj^,*,l)  a submonoid.  Furthermore, 
each  (Zj^,  + ,0)  is  an  abelian  group,  or  a group  for  which  + 
is  commutative  as  well  as  associative  in  Zj^,  and  is 
referred  to  as  the  additive  group  of  the  ring  ( Zj^ , + , ’ , 0 , 1 ) . 
Each  triple  (Zj^,’,l)  is  called  the  multiplicative  monoid  of 
the  ring. 

A ring  ( S , <|)j^ , <i)2  , U2  ) is  a structure  comprising 

a nonvacuous,  S,  set  together  with  two  compositions  '♦>2 

in  S and  two  distinguished  unit  elements  Uj^,U2sS  such  that 
( S , <|)j^  ,Uj^ ) forms  an  abelian  group  (a  group  for  which  <}>j^  is 
commutative  as  well  as  associative  in  S)  and  (S,<j>2,U2) 
forms  a monoid  and  for  all  a,b,ceS,  a<|>2  ( t»<l'j^c ) =a<t>2b<j>j^a<j>2C 
and  ( b(|>j^c ) <|)2=b<|)2a<l>2c4>23  • example,  the  integers  form  a 

ring  over  addition  and  subtraction  called  the  Gaussian  ring 
of  integers  and  denoted  (Z,+,’,0,1).  If  the  nonvacuous 
elements  form  an  abelian  subgroup  of  the  monoid  (S,<t>2,U2), 
then  ( S , <t>j^ , <|>2  f ,U2  ) defines  a field . More  precisely 
then,  the  residue  sets  ( Zj^ , + , ' , 0 , 1 ) , k=l,2,...,L  each  form 
a field. 

If  S^,  S2,--.,Sj^  each  form  respective  monoids  (groups) 
under  <|>2  / • • • / / then  their  direct  product  S^XS2X . . . 

also  forms  a monoid  (group)  whose  elements  are  the  set  of 
L-tuples,  each  denoted  (a^,  a2,..-,aL)  a^eZ^ 
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i=l,2,...,L  under  the  following  composition: 


/ 


(A-3) 


An  important  RNS  property  follows  from  this  relation.  That 
is,  if  a,b,ceZj^  and  c=a<j>b  where  a->(  a^  , a2  , . . • , ) and 
b->(  bj^ , b2  , . . . , bj^ ) correspond  to  the  L-tuple  of  residues  each 
taken  modulo  pj^,  k = l,2,...,L,  then  c-»(  c^  , C2  , . . . , Cj^ ) where 
C]^=(aj^<|)bj^)modp]^,  k=l,2,...,L. 

In  words,  addition  and  multiplication  can  be  computed 
via  L concurrent  modulo  Pj^  operations  without  the  need  for 
interdigit  communication.  As  a result,  the  RNS  is  referred 
to  as  a carry-free  algebraic  system  where  each  operation 
aj^ijibj^  is  implemented  as  a table-lookup  using  high-speed 
semiconductor  memory. 

As  a final  point,  in  the  RNS  the  dynamic  range  (given 
by  A-2)  may  be  exceeded  during  processing  as  long  as  the 
final  output  does  not  exceed  the  prespecified  range.  To 
illustrate,  suppose  the  operation  xy-z  is  to  be  performed 
on  the  following  complex  data:  x=8+j3,  y=3+j2,  and 
z=10+jll.  Direct  multiplication  of  x and  y yields  the 
complex  product  ( 8+ j 3 ) ( 3+ j 2 ) =18+ j 25  from  which  10+jll  is 
subtracted  to  produce  a final  value  of  8+jl4.  If  taken 
congruent  modulo  p where  p=17,  the  product  is 
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<18+ j25>j^^sl  + j8  which  in  a single  modulus  system  (L=l) 
would  produce  ambiguous  output  as  a stand-alone  task. 
Coupled  with  the  operation  -zmodulop  (z=10+jll),  however, 
the  output  <-9-jl3>^.ys8+ jl4  is  obtained,  which  is  the 
expected  output. 


APPENDIX  B 

OPTICAL  TRANSFORMATION  PRINCIPLES 


The  transformation  properties  thin  of  lenses  are 
well-known  [74]  and  are  given  briefly  below.  Let  a^(x,z) 
denote  the  input  optical  wave  amplitude  of  a time  harmonic 
wave  travelling  in  the  z direction  (and  independent  of  y) 
and  incident  on  a thin  lens.  A standard  assumption  of  a 
thin  lens  is  that  the  transverse  displacement  of  an  optical 
beam  passing  through  the  lens  be  small  such  that  a ray 
incident  on  the  front  surface  of  the  lens  at  x exits  from 
the  back  surface  of  the  lens  at  x with  a phase  delay  of  the 
form 


(|)(x)=knA(x)-k[A(0)-A(x)  ] 


(B-1) 


where  n is  the  refractive  index  of  the  lens  material,  A(x) 
denotes  lens  thickness,  and  k=2n/X  where  X is  the 
wavelength  of  the  optical  wave.  In  this  case,  the  complex 
amplitude  of  the  emergent  optical  wave  is  a phase  delayed 


172 


173 


version  of  the  incident  wave  amplitude  which  satisfies  the 
following  relation 


(x)=e“^  a,  (x)l(x) 


(B-2) 


(1  inside  the  lens  aperature 
0 elsewhere 


for  k=2n/X,  where  X is  the  free  space  wavelength  of  the 
incident  field,  n the  refractive  index  of  the  lens 
material,  and  f the  lens  focal  length  as  a function  of  n 
and  lens  curvature.  (For  a thin  lens  inclusion  of  the  lens 
thickness  imparts  only  a constant  phase  term  of  the  form 
g-]nA(0)  where  A(0)  denotes  the  maximum  thickness  of  the 
lens  and  is  typically  small.) 

The  Fresnel  propagation  approximation  can  be  applied 

at  the  back  focal  plane  of  the  lens.  The  Fresnel 

propagation  approximation  is  essentially  a first  order 

binomial  approximation  on  Huygen's  principle,  which  states 

that  each  point  on  a wavefront  may  be  reguarded  as  a new 

source  of  waves.  Each  wave  element  dx  contributes  to  the 

complex  amplitude  of  a wavefront  at  point  P the  equivalent 

n k R 

of  ( EQa^  ( X )/XR)  e-*  dx  producing  a complex  amplitude 

distribution  which  expresses  as 
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OD 

a(Xp)  = (EQ/X)J(l/R)ag(x)e^’^^dx 
— 00 


(B-3.a) 


where 


R=[f^+(x-Xp)^]^/^=f [l+( (x-Xp)/X)2]^/2 
= Zp(l+(l/2)  { (x-Xp)/Zp)^-(l/8)  ( (x-Xp)/Zp)'^+.  . . ) (B-3.b) 


for  Eq  a constant  and  R assumed  >>  X.  At  the  back  focal 
plane  of  the  transform  lens,  a first  order  approximation 
can  be  applied  to  R in  the  exponent  where  Zp=f  • In  the 
denominator  small  changes  in  R are  not  as  important  as  in 
the  exponent  and  R can  be  approximated  by  f.  Equation 
(B-3)  reduces  to 

CO 

_co 

where 

t(Xj).(Ej,/jX£)e3“e3''='f/2f 

for  Xj/Xf  termed  the  spatial  frequency  coordinate.  The 
term  e-*  is  a constant  phase  term  which  can  be  neglected 
as  it  equates  to  unity  when  the  square  modulus  is  taken  to 
compute  light  intensity.  Substituting  (B-2)  in  for  3l^{x) 
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yields 

00 

a(x£)=t(x^)Ja.  (x)l(x)e"^^'^''^''f/^^^dx.  (B-5) 

— 00 

If  a^(x)  denotes  a time  harmonic  transmittance  wave 
amplitude  incident  on  the  lens  imposed  by  a diffraction 
grating  s(x)  on  a uniform  monochromatic  light  source  a 
distance  dO  in  front  of  the  lens,  the  following  relation 
holds [ 1 ] 

(B-6) 

where  S(Xj)  is  the  Fourier  spectrum  of  s(x)  and  A^(x^)  is 
the  Fourier  spectrum  of  a^(x).  Using  this  relation  (B— 5) 
becomes 


a(x£)  = (EQ/Xf  )t' (Xj)  Js(x)l(x)e  j 2 ixx  ( x ^/Xf  ) 

.00 

where 

For  dg=f,  t'(xj)=l  and  (B-7)  is  the  Fourier  transform  in 
x^/Xf  of  the  diffraction  function,  s(x),  subtended  by  the 
lens  aperture. 


APPENDIX  C 
GLOSSARY 


AF  (ambiguity  function) 

A two-dimensional  representation  of  signal  delay  and 
doppler  shift  common  to  radar  and  sonar  applications. 

BT  (bandwidth-time  product) 

The  product  of  the  effective  frequency  bandwidth,  the 
frequency  band  over  which  the  amplitude  of  the  normalized 
spectrum  exceeds  exp[-n/4],  and  the  effective  time 
duration,  the  time  duration  over  which  the  amplitude  of  the 
normalized  signal  exceeds  exp[-n/4]. 

P (dispersion  index) 

An  indicator  of  the  spectral  concentration  of  a 
signal.  For  FM  signals,  the  dispersion  index  is  defined  as 
the  product  of  the  period  of  the  modulating  signal  and  its 
amplitude  where  the  amplitude  is  the  maximum  frequency 
deviation  from  the  carrier  frequency  of  the  FM  signal. 
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CRNS  (complex  residue  numbering  system) 

A finite  ring  of  complex  integers  where  each  element 
of  the  ring  is  uniquely  represented  by  a k-tuple  of  complex 
residues  taken  modulo  a set  of  relative  primes. 

DFT  (discrete  Fourier  transform) 

The  N samples  taken  at  w=2nk/N  of  the  decompostion  of 
an  N point  data  sequence  into  linear  combinations  of 

DWD  (discrete  Wigner  distribution) 

The  discrete  form  of  the  Wigner  distribution  in  both 
time  and  frequency  and  used  in  digital  signal  processing 
applications . 

FFT  (fast  Fourier  transform) 


Any  member  of 

a 

class  of 

algorithms  used 

to 

efficiently 

compute 

the 

discrete 

Fourier 

transform 

by 

decomposing 

either 

the 

input  or 

output 

sequence 

into 

integer  powers. 

QRNS  (quadratic  residue  number  system) 

A complex  residue  representation  where  each  prime,  p, 

2 

is  chosen  such  that  the  quadratic  i s-lraodulop  has  a 
solution  which  allows  an  isomorphism  from  each  residue 
representaion  to  a two-tuple  representation  generated  by 
the  product  set  of  integers  modulo  p with  itself. 
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RNS  (residue  number  system) 

A finite  ring  of  integers  where  each  element  of  the 
ring  is  uniquely  represented  by  a k-tuple  of  residues  taken 
modulo  a set  of  relative  primes. 

SM-QRNS  (single-modulus  quadratic  residue  number 
system) 

The  quadratic  residue  number  system  for  the  case  of  a 
single  modulus  (i.e.  one  prime). 

spDWD  ( smoothed-pseudo  discrete  Wigner  Distribution) 

The  discrete  form  of  the  weighted  Wigner  distribution 
in  both  time  and  frequency  and  used  in  digital  signal 
processing  applications. 

STFT  (short-time  Fourier  transform) 

The  Fourier  transform  of  a windowed  signal  where 
stationarity  is  postulated  over  the  length  of  the  window  as 
the  window  slides  along  the  time  axis. 

STP  (short-time  periodogram) 

The  square  magnitude  of  the  short-time  Fourier 


transform. 
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WD  (Wigner  distribution) 

A two-dimensional  representation  of  signal 
concentration  over  time  and  frequency  which  maintains  the 
characteristic  signal  properties  of  the  one-dimensional 
time  and  frequency  domains  yet  attains  negative  values. 

wWD  (weighted  Wigner  distribution) 

The  resulting  representation  after  moving  averages,  in 
the  form  of  windowing  and  filtering  operations,  are  applied 
to  the  Wigner  distribution. 


APPENDIX  D 
PROGRAM  LISTING 


C COS 

generates  sinusoidal  data,  then  makes 

compatible 

INTEGER*4  ICOMMON(128) 

INTEGER*2  ITRY(256) 

REAL  F0,TPI,RNU 
INTEGER  M 

CHARACTER*8  NAMEI,NAMEO 
PRINT  *,  'enter  FO' 

READ  *,  FO 
WRITE(6,100) 

100  FORMATdx,  'enter  the  sampled  data  file') 

READ( 5,110 )NAMEI 
110  FORMAT d8) 

WRITE(6,120) 

120  FORMATdx,  'enter  the  cos  file  name') 

READ{ 5 , 130 )NAMEO 
130  FORMAT (A8) 


it  ILS 


OPEN ( UNIT=1 , NAME=NAMEI , STATUS= ' OLD ' , FORM= ' UNFORMATTED ' 

1 READONLY)  ' 


OPEN ( UNI T=2 , NAME=NAMEO , STATUS= ' NEW ' , FORM= ' UNFORMATTED ' 

1 ACCESS='DIRECT' ,RECL=128 ) 

C initialize  parameters 

period  (eg  1536.)' 

i\CiAU  ^ f RNU 
TPI=8 . *ATAN(1 . ) 

IREC=1 


M=-256 

C set  header  block  for  ils  compatibility 
READd  ) (ICOMMONd  ) , 1 = 1,128) 

WRITE( 2, REC=IREC) (ICOMMONd ) , 1=1,128) 
250  READd  , END=900  ) ( ITRY(  I ) , 1 = 1,256) 

C Generate  cos  data 
M=M+256 


DO  200  J=l,256 

ITRY( J)=1000 . *COS(TPI*F0*(M+J)/RNU) 
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200 

900 


CONTINUE 

IREC=IREC+1 

WRITE(2,REC=IREC) ( ITRY( I ) , I 
GO  TO  250 
CONTINUE 
CLOSE( 1 ) 

CLOSE( 2 ) 

STOP 

END 


=1,256) 
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C DWD 

C This  program  takes  ILS  compatible  sampled  data,  generates 
the  Wigner 

C inner  product,  then  outputs  the  IP  as  an  ILS  compatible 
data  file; 

C in  this  form,  when  passed  through  the  ILS  program  SDI  the 
output  is 

C the  DWD  of  the  original  data. 

INTEGER*4  ICOMMON(128) 

INTEGER*2  ITRY(256) 

REAL  FITRY(256) 

CHARACTER*8  NAMEI, NAMED 
WRITE(6,100) 

100  FORMATdX,  'enter  the  sampled  data  file') 

READ( 5 , 110 )NAMEI 
110  FORMAT(A8) 

WRITE(6,120) 

120  FORMATdX,  'enter  the  DWD  IP  file  name') 

READ( 5,130)NAMEO 
130  FORMAT (A8) 

OPEN ( UNIT=1 , NAME=NAMEI , STATUS= ' OLD ' , FORM= ' UNFORMATTED ' 

1 READONLY) 

OPEN ( UNIT=2 , NAME=NAMEO , STATUS= ' NEW ' , FORM= ' UNFORMATTED ' 

1 ACCESS='DIRECT' ,RECL=128 ) 

C initialize  parameters 
IREC=1 

C set  header  block  for  ils  compatibility 
READd  ) (ICOMMONd  ) , 1 = 1,128) 

WRITE( 2,REC=IREC) (ICOMMONd ) , 1=1,128) 

250  READ( 1 , END=900 ) (ITRYd ) , 1=1,256) 

C generate  DWD  inner  product 
DO  200  J=l,128 

FITRY( J)=FLOAT(ITRY( J) )*FLOAT( ITRY(257-J) )/1000. 

FITRY( 257-J)=FITRY( J) 

200  CONTINUE 

DO  300  J=l,256 

ITRY( J)=NINT( FITRY( J) ) 

300  CONTINUE 

IREC=IREC+1 

WRITE( 2 ,REC=IREC) ( ITRY( I ) , 1=1,256) 

GO  TO  250 
900  CONTINUE 

CLOSEd  ) 

CLOSE( 2 ) 

STOP 

END 
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C SQRMOD 

C This  program  generates  data  from  square  pulse  modulation 
C on  a carrier,  then  makes  it  ILS  compatible 
INTEGER*4  ICOMMON(128) 

INTEGER*2  ITRY(256) 

REAL  F0,TPI,RNU 
INTEGER  M 

CHARACTER* 8 NAME I, NAMED 
WRITE(6,100) 

100  FORMAT(ix,  'enter  the  sampled  data  file') 

READ( 5 , 110 )NAMEI 
110  FORMAT (A8) 

WRITE(6,120) 

120  FORMAT(lX,  'enter  the  square  file  name') 

READ( 5 , 130 ) NAMED 
130  FORMAT (A8) 

OPEN( UNIT=1 , NAME=NAMEI , STATUS= ' OLD ' , FORM= ' UNFORMATTED' , 

1 READONLY) 

OPEN( UNIT=2 ,NAME=NAMEO , STATUS= ' NEW' , FORM= ' UNFORMATTED ' , 

1 ACCESS='DIRECT' ,RECL=128 ) 

C initialize  parameters 

PRINT  *,  'enter  the  carrier  frequency' 

READ  *,  FO 

PRINT  *,  'enter  data  points  per  period  (eq 

1536.)'  V y 

READ  *,RNU 
TPI=8 . *ATAN( 1 . ) 

IREC=1 

M=-256 

C set  header  block  for  ils  compatibility 
READ( 1 ) ( ICOMMON( I ) , 1=1,128) 

WRITE( 2 ,REC=IREC) ( I COMMON ( I ) , 1=1,128) 

250  READ( 1 , END=900 ) ( ITRY( I ) , 1=1,256) 

C Generate  square  pulse  data 
M=M+256 

DO  200  J=l,256 
VAL=M+J 
HNU=RNU/2 

983  IF(VAL  .GT.  RNU)THEN 

VAL=VAL-RNU 
GO  TO  983 
ELSE 

IF(VAL  .LE.  HNU)THEN 
ITRY( J ) =1000 . *COS ( TPI*F0* ( M+J ) /RNU ) 

ELSE 

ITRY( J)=0 . 

END  IF 
END  IF 

200  CONTINUE 

IREC=IREC+1 
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WRITE(2,REC=IREC) ( ITRY( I) , I 
GO  TO  250 
CONTINUE 
CLOSE( 1 ) 

CLOSE( 2 ) 

STOP 

END 


=1,256) 
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C DWD64 

C This  program  takes  ILS  compatible  sampled  data,  generates 
the  Wigner 

C inner  product,  then  outputs  the  IP  as  an  ILS  compatible 
data  file; 

C in  this  form,  when  passed  through  the  ILS  program  SDI  the 
output  is 

C the  DWD  of  the  original  data. 

INTEGER*4  ICOMMON(128) 

INTEGER*2  ITRY(256) 

REAL  FITRY(256) 

CHARACTER*8  NAMEI, NAMED 
WRITE( 6 , 100 ) 

100  FORMAT(lX,  'enter  the  sampled  data  file') 

READ( 5 , 110 )NAMEI 
110  FORMAT (A8) 

WRITE(6,120) 

120  FORMATdX,  'enter  the  inner  product  file  name') 

READ( 5 , 130 )NAMEO 
130  FORMAT (A8) 

OPEN ( UNIT=1 , NAME=NAMEI , STATUS= ' OLD ' , FORM= ' UNFORMATTED' , 

1 READONLY) 

OPEN( UNIT=2 , NAME=NAMEO , STATUS= ' NEW' , FORM= ' UNFORMATTED ' , 

1 ACCESS='DIRECT' ,RECL=128 ) 

C initialize  parameters 
IREC=1 

C set  header  block  for  ils  compatibility 
READ( 1 ) ( ICOMMON( I ) , 1=1,128) 

WRITE( 2 ,REC=IREC) ( ICOMMON( I ) , 1=1 , 128 ) 

250  READ( 1 , END=900 ) ( ITRY( I ) , 1=1,256) 

C generate  DWD  inner  product 
DO  200  J=l,32 

FITRY( J )=FLOAT( ITRY( J ) ) * FLOAT ( I TRY ( 65- J ) )/1000 . 

FITRY( J+64 )=FLOAT( ITRY( J+64) ) *FLOAT( ITRY( 129-J ) )/1000, 

FITRY( J+128 )=FLOAT( ITRY( J+128 ) ) *FLOAT( ITRY( 193-J ) )/1000 . 

FITRY( J+192 )=FLOAT( ITRY( J+192 ) )* FLOAT ( ITRY( 257-J ) )/1000 . 
FITRY( 65-J)=FITRY( J) 

FITRY( 129-J )=FI TRY ( J+64 ) 

FITRY( 193-J)=FITRY( J+128 ) 

FITRY( 257-J )=FITRY( J+192 ) 

200  CONTINUE 

DO  300  J=l,256 
ITRY( J )=NINT( FITRY( J ) ) 

300  CONTINUE 

IREC=IREC+1 

WRITE( 2,REC=IREC) ( ITRY( I ) , 1=1,256) 

GO  TO  250 
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900 


CONTINUE 
CLOSE( 1 ) 
CLOSE( 2 ) 
STOP 
END 
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C DWDW 

C This  program  takes  ILS  compatible  sampled  data, 

C windows  it,  generates  the  Wigner 

C inner  product,  then  outputs  the  IP  as  an  ILS  compatible 
data  file; 

C in  this  form,  when  passed  through  the  ILS  program  SDI  the 
output  is 

C the  DWD  of  the  original  data. 

INTEGER*4  ICOMMON(128) 

INTEGER*2  ITRY(256) 

REAL  FITRY(256),  WNDW( 128 ) , PITRY( 256 ) 

CHARACTER*8  NAMEI,NAMEO 
WRITE( 6 ,100 ) 

100  FORMAT(lX,  'enter  the  sampled  data  file') 

READ( 5 , 110 )NAMEI 
110  FORMAT (A8) 

WRITE(6,120) 

120  FORMATdX,  'enter  the  DWD  IP  file  name') 

READ( 5 , 130 )NAMEO 
130  FORMAT (A8) 

OPEN ( UNI T=1 , NAME=NAMEI , STATUS= ' OLD ' , FORM= ' UNFORMATTED ' , 

1 READONLY) 

OPEN(UNIT=2,NAME=NAMEO,STATUS='NEW' , FORM= ' UNFORMATTED ' , 

1 ACCESS='DIRECT' ,RECL=128 ) 

C initialize  parameters 
IREC=1 

READ( 3, * ) (WNDW( J) , J=l,128) 

C set  header  block  for  ils  compatibility 
READ( 1 ) ( ICOMMON( I ) , 1=1,128) 

WRITE( 2,REC=IREC) (ICOMMON(I) , 1=1,128) 

250  READ(1,END=900) (ITRY(I) , 1=1,256) 

C generate  DWD  inner  product 
DO  200  J=l,128 

PITRY( J)=FLOAT(ITRY( J) )*WNDW( J) 
PITRY(257-J)=FLOAT( ITRY(257-J) )*WNDW( J) 

FITRY( J)=PITRY( J) *PITRY( 257-J)/1000 . 
FITRY(257-J)=FITRY( J) 

200  CONTINUE 

DO  300  J=l,256 
ITRY( J)=NINT(FITRY( J) ) 

300  CONTINUE 

IREC=IREC+1 

WRITE( 2 ,REC=IREC ) ( ITRY( I ) , 1=1,256) 

GO  TO  250 
900  CONTINUE 

CLOSE( 1 ) 

CLOSE( 2 ) 

CLOSE( 3 ) 

STOP 

END 
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C DWD64 

C This  program  takes  ILS  compatible  sampled  data,  generates 
the  Wigner 

C inner  product,  then  outputs  the  IP  as  an  ILS  compatible 
data  file; 

C in  this  form,  when  passed  through  the  ILS  program  SDI  the 
output  is 

C the  DWD  of  the  original  data. 

INTEGER*4  ICOMMON(128) 

INTEGER*2  ITRY(256) 

REAL  FITRY(256) ,PITRY(256) ,WNDW(32) 

CHARACTER*8  NAMEI, NAMED 
WRITE(6,100) 

100  FORMATdX,  'enter  the  sampled  data  file') 

READ( 5, 110) NAMEI 
110  FORMAT (A8) 

WRITE(6,120) 

120  FORMAT(lX,  'enter  the  inner  product  file  name') 

READ (5,130) NAMED 
130  FORMAT (A8) 

OPEN( UNIT=1 ,NAME=NAMEI , STATUS='OLD' , FORM= ' UNFORMATTED' 

1 READONLY) 

OPEN(UNIT=2,NAME=NAMEO,STATUS='NEW' , FORM= ' UNFORMATTED ' , 

1 ACCESS='DIRECT' ,RECL=128 ) 

C initialize  parameters 
IREC=1 

READ(3,*) (WNDW( J) , J=l,32) 

C set  header  block  for  ils  compatibility 
READ( 1 ) ( ICOMMON( I ) , 1=1,128) 

WRITE( 2 ,REC=IREC) ( I COMMON ( I ) , 1=1,128) 

250  READ( 1 ,END=900 ) ( ITRY( I ) , 1=1,256) 

C generate  DWD  inner  product 
DO  200  J=l,32 

PITRY( J)=FLOAT( ITRY( J ) ) *WNDW( J ) 

PITRY( J+64 )=FLOAT( ITRY( J+64 ) ) *WNDW( J ) 

PITRY( J+128 )=FLOAT( ITRY( J+128 ) ) *WNDW( J) 

PITRY( J+192)=FLOAT( ITRY( J+192) ) *WNDW( J) 

PITRY( 65-J )=FLOAT( ITRY( 65-J ) ) *WNDW( J ) 
PITRY(129-J)=FLOAT( ITRY( 129-J) ) *WNDW( J) 
PITRY(193-J)=FLOAT( ITRY(193-J) )*WNDW( J) 
PITRY(257-J)=FLOAT( ITRY(257-J) )*WNDW( J) 

FITRY( J )=PITRY( J) *PITRY( 65-J )/1000 . 

FITRY( J+64 )=PITRY( J+64 ) *PITRY( 129-J )/1000 . 

FITRY( J+128 )=PITRY( J+128 ) *PITRY( 193-J )/1000 . 

FITRY( J+192 )=PI TRY ( J+192 ) *PITRY( 257-J )/1000 . 

FITRY( 65-J)=FITRY( J ) 

FITRY( 1 29-J)=FI TRY ( J+64 ) 

FI TRY ( 193-J )=FITRY( J+128) 

FITRY( 257-J )=FI TRY ( J+192 ) 

200  CONTINUE 
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300 

900 


DO  300  J=l,256 
ITRY( J)=NINT( FITRY(  J) ) 
CONTINUE 
IREC=IREC+1 

WRITE( 2 ,REC=IREC) ( ITRY(  I ) , I 
GO  TO  250 
CONTINUE 
CLOSE( 1 ) 

CLOSE( 2 ) 

STOP 

END 


=1,256) 
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C SCIRMOD2 

C This  program  generates  data  from  sin  pulse  modulation 
C on  a carrier,  then  makes  it  ILS  compatible 
INTEGER*4  ICOMMON(128) 

INTEGER*2  ITRY(256) 

REAL  FO,TPI,RNU,FMOD 
INTEGER  M 

CHARACTER*8  NAMEI, NAMED 
WRITE(6,100) 

100  FORMATdX,  'enter  the  sampled  data  file') 

READ( 5,110)NAMEI 
110  FORMAT(A8) 

WRITE(6,120) 

120  FORMAT(lX,  'enter  the  sine  pulse  file  name') 

READ( 5, 130 ) NAMED 
130  FDRMAT(A8) 

DPEN ( UNIT=1 , NAME=NAMEI , STATUS= ' DLD ' , FDRM= ' UNFDRMATTED ' 

1 READDNLY) 


DPEN ( UNI T=2 , NAME=NAMED , STATUS= ' NEW ' , FDRM= ' UNFDRMATTED ' 
_1  ACCESS='DIRECT' ,RECL=128 ) 

C initialize  parameters 


1536.  ) ' 


PRINT  *,  'enter  the  carrier  frequency' 

READ  *,  FO 

PRINT  *,  'enter  data  points  per  period  (eg 


READ  * , RNU 
TPI=8 . *ATAN( 1 . ) 
PI=TPI/2 
IREC=1 


M=-256 

C set  header  block  for  ils  compatibility 
READ( 1 ) ( ICDMMDN( I ) , 1=1,128) 

WRITE ( 2 ,REC=IREC) ( ICDMMDN( I ) , 1=1,128) 
250  READ( 1 ,END=900 ) ( ITRY( I ) , 1=1,256) 

C Generate  sine  pulse  data 
M=M+256 


t 


DD  200  J=l,256 
VAL=M+J 


HNU=RNU/2 

983  IF(VAL  .GT.  RNU) THEN 

VAL=VAL-RNU 
GD  TD  983 
ELSE 

FMDD=1000 . *CDS( PI*VAL/RNU) 

ITRY( J)=FMDD*CDS(TPI*F0*(M+J) /RNU) 
END  IF 

200  CDNTINUE 

IREC=IREC+1 

WRITE( 2 ,REC=IREC) ( ITRY( I ) , 1=1,256) 

GD  TD  250 
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900 


CONTINUE 
CLOSE( 1 ) 
CLOSE( 2 ) 
STOP 
END 
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C FMCOS 

C This  program  generates  data  from  FM  sin  pulse  modulation 
C on  a carrier,  then  makes  it  ILS  compatible 
INTEGER*4  ICOMMON(128) 

INTEGER*2  ITRY(256) 

REAL  FO ,TPI ,RNU, FMOD 
INTEGER  n 

CHARACTER*8  NAMEI, NAMED 
WRITE(6,100) 

100  FORMAT(ix,  'enter  the  sampled  data  file') 

READ( 5, 110 )NAMEI 
110  FORMAT (A8) 

WRITE(6,120) 

120  FORMATdX,  'enter  the  FM  file  name') 

READ( 5 , 130 ) NAMED 
130  FORMAT (A8) 

OPEN ( UNIT=1 , NAME=NAMEI , STATUS= ' OLD ' , FORM= ' UNFORMATTED ' , 

1 READONLY) 

OPEN(UNIT=2,NAME=NAMEO,STATUS='NEW' , FORM= ' UNFORMATTED ' , 

1 ACCESS='DIRECT' ,RECL=128 ) 

C initialize  parameters 

PRINT  *,  'enter  the  carrier  frequency' 

READ  *,  FO 

PRINT  *,  'enter  data  points  per  period  (eq 
1536  .)'  V y 

READ  * , RNU 

PRINT  *,  'enter  GAIN  value' 

READ  *,  GAIN 
TPI=8 . *ATAN( 1 . ) 

PI=TPI/2 

IREC=1 

M=-256 

C set  header  block  for  ils  compatibility 
READ( 1 ) ( ICOMMON( I ) , 1=1,128) 

WRITE( 2 ,REC=IREC ) ( I COMMON ( I ) , 1=1,128) 

250  READ( 1 , END=900 ) ( ITRY( I ) , 1=1,256) 

C Generate  FM  data 
M=M+256 

DO  200  J=l,256 
VAL=M+J 

983  IF(VAL  .GT.  RNU) THEN 

VAL=VAL-RNU 
GO  TO  983 
ELSE 

FMOD=l .+COS( PI*VAL/RNU) 

ITRY( J)=GAIN*COS( TPI*F0*FMOD*VAL/RNU) 

END  IF 

200  CONTINUE 

IREC=IREC+1 

WRITE( 2,REC=IREC) ( ITRY( I ) , 1=1,256) 
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900 


GO  TO  250 
CONTINUE 
CLOSE( 1 ) 
CLOSE( 2 ) 
STOP 
END 


C CHRP 

C This  program  generates  chirp  data,  then  makes  it  ILS 
compatible 

INTEGER*4  ICOMMON(128) 

INTEGER*2  ITRY(256) 

REAL  F0,TPI,RNU 
INTEGER  M 

CHARACTER*8  NAMEI, NAMED 
PRINT  *,  'enter  FO' 

READ  * , FO 
WRITE(6,100) 

100  FORMAT(lX,  'enter  the  sampled  data  file') 

READ( 5 , 110 ) NAMEI 
110  FORMAT (A8) 

WRITE(6,120) 

120  FORMAT(lX,  'enter  the  chirp  file  name') 

READ( 5, 130) NAMED 
130  FORMAT (A8) 

OPEN(UNIT=l ,NAME=NAMEI , STATUS='OLD' , FORM= ' UNFORMATTED ' 
1 READONLY) 

OPEN ( UNI T=2 , NAME=NAMEO , STATUS= ' NEW' , FORM= ' UNFORMATTED' 
1 ACCESS='DIRECT' ,RECL=128 ) 

C initialize  parameters 

PRINT  *,  'enter  data  points  per  period  (eq 
1536.)'  V y 

READ  * , RNU 
TPI=8 . *ATAN( 1 . ) 

IREC=1 

M=-256 

C set  header  block  for  ils  compatibility 
READ( 1 ) ( ICOMMON( I ) , 1=1,128) 

WRITE( 2 ,REC=IREC ) ( I COMMON ( I ) , 1=1,128) 

250  READ( 1 , END=900 ) ( ITRY( I ) , 1=1,256) 

C Generate  chirp  data 
M=M+256 

DO  200  J=l,256 

ITRY( J)=1000.*COS(TPI*F0*( ( M+ J ) /RNU ) * * 2/2 . ) 
200  CONTINUE 

IREC=IREC+1 

WRITE( 2 ,REC=IREC) ( ITRY( I ) , 1=1,256) 

GO  TO  250 
900  CONTINUE 

CLOSE( 1 ) 

CLOSE( 2 ) 

STOP 

END 


C FMPULS 

C This  program  generates  data  from  square  pulse  FM 
modulation 

C on  a carrier,  then  makes  it  ILS  compatible 
INTEGER*4  ICOMMON(128) 

INTEGER*2  ITRY(256) 

REAL  FO ,TPI ,RNU, FMOD 
INTEGER  M 

CHARACTER*8  NAMEI, NAMED 
WRITE(6,100) 

100  FORMAT(lX,  'enter  the  sampled  data  file') 

READ( 5, 110 )NAMEI 
110  FORMAT (A8) 

WRITE( 6 , 120 ) 

120  FORMAT(lX,  'enter  the  FM  file  name') 

READ( 5 , 130 )NAMEO 
130  FORMAT (A8) 

OPEN ( UNIT=1 , NAME=NAMEI , STATUS= ' OLD ' , FORM= ' UNFORMATTED 
1 READONLY) 

OPEN ( UNIT=2 , NAME=NAMEO , STATUS= ' NEW' , FORM= ' UNFORMATTED 
1 ACCESS='DIRECT' ,RECL=128 ) 

C initialize  parameters 

PRINT  *,  'enter  the  carrier  frequency' 

READ  *,  FO 

PRINT  *,  'enter  data  points  per  period  (eq 

1536.)'  V y 

READ  *,RNU 

PRINT  *,  'enter  GAIN' 

READ  *,  GAIN 
TPI=8 . *ATAN( 1 . ) 

PI=TPI/2 

IREC=1 

M=-256 

C set  header  block  for  ils  compatibility 
READ( 1 ) ( ICOMMON( I ) , 1=1,128) 

WRITE( 2 ,REC=IREC ) ( I COMMON ( I ) , 1=1,128) 

250  READ( 1 ,END=900 ) ( ITRY( I ) , 1=1,256) 

C Generate  pulse  FM  data 
M=M+256 

DO  200  J=l,256 
VAL=M+J 
HNU=RNU/2 

983  IF(VAL  .GT.  RNU)THEN 

VAL=VAL-RNU 
GO  TO  983 
ELSE 

IF(VAL  .LE.  HNU)THEN 
FMOD=1000 . 

ELSE 

FMOD=500 . 
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200 

900 


END  IF 
END  IF 


I TRY ( J ) =GAIN*COS ( TPI *F0 *FMOD*VAL/RNU ) 
CONTINUE 


IREC=IREC+1 

WRITE( 2,REC=IREC) ( ITRY( I ) , 1=1,256 ) 
GO  TO  250 
CONTINUE 
CLOSE( 1 ) 

CLOSE( 2 ) 

STOP 

END 
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C cos2 

generates  data  from  the  sum  of  two 

o^nusoxos  f 

C then  makes  it  ILS  compatible 
INTEGER*4  ICOMMON(128) 

INTEGER*2  ITRY(256) 

REAL  F0,TPI,RNU,F1 
INTEGER  M 

CHARACTER* 8 NAMEI, NAMED 
PRINT  *,  'enter  FO' 

READ  *,  FO 

PRINT  *,  'enter  Fl' 

READ  *,  Fl 

WRITE(6,100) 

FORMATdX,  'enter  the  sampled  data  file') 
READ(5,110)NAMEI  ’ 

FORMAT (A8 ) 

WRITE(6,120) 

FORMAT (A8 ) 


100 

110 

120 

130 


OPEN( UNIT-I^NAME-NAMEI , STATUS- ' OLD ■ , FORM- ’ UNFORMATTED ’ , 

0PEN(UNIT=2  NAME-NAHEO.STATUS.'NEW' .FORM-'UNFORMATTED' , 

_ 1 ACCESS='DIRECT' ,RECL=128 ) 

C initialize  parameters 

1536  .)'  *'  'enter  data  points  per  period  (eg 

READ  *,RNU 
TPI=8 . *ATAN( 1 . ) 

IREC=1 

M=-256 

C set  header  block  for  ils  compatibility 
READd  ) ( ICOMMON(  I ) , 1 = 1,128) 

WRITE(2,REC=IREC)  (ICOMMONd)  , 1 = 1,128) 

250  READd , END=900 ) (ITRYCI ) , 1=1,256) 

C Generate  cos2  data 
M=M+256 

DO  200  J=l,256 

ITRY(  J)  = 1000  .(  COS(TPI»FOMM+J)/RNU)+COS(TPI*F1MM+J)/RNU)  ) 


900 


CONTINUE 
IREC=IREC+1 

WRITE(2,REC=IREC)  (ITRYd)  , 1 = 1,2  56) 
GO  TO  250  > 

CONTINUE 
CLOSEd  ) 

CLOSE( 2 ) 

STOP 

END 
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C SPDWD 

C This  program  takes  ILS  compatible  sampled  data,  generates 
the  SPDWD 

j P>^oduct,  then  outputs  the  IP  as  an  ILS  compatible 

data  file; 

C in  this  form,  when  passed  through  the  ILS  program  SDI  the 
output  is 


C the 


100 

110 

120 

130 


SPDWD  of  the  original  data. 

INTEGER*4  ICOMMON(128) 

INTEGER*2  ITRY(256) ,IOUT(256) , SMV 
REAL  TTRY( 768 ) ,RIP 
CHARACTER*8  NAMEI,NAMEO 
WRITE(6,100) 

FORMATdX,  'enter  the  sampled  data  file') 
READ( 5,110 ) NAME I 
FORMAT (A8 ) 

WRITE(6,120) 

FORMATdX,  'enter  the  SPDWD  IP  file  name') 
READ( 5 , 130 )NAMEO 
FORMAT (A8 ) 


OPEN(UNIT=l ,NAME=NAMEI , STATUS='OLD' , FORM= ' UNFORMATTED ' 

1 READONLY)  ' 


OPEN ( UNIT=2 , NAME=NAMEO , STATUS= ' NEW' , FORM= ' UNFORMATTED ' 

_1_  ACCESS='DIRECT' ,RECL=128 ) ' 

C initialize  parameters 
IREC=1 
WRITE(6,11) 

FORMAT(lX,  'enter  the  smoothing  parameter  value') 

READ(5,*)SMV 

RSMV=REAL( SMV) 

WRITE(6,12) 

12  FORMAT(lX,  'enter  the  necessary  magscale  for  ILS 

compatibility' ) 

READ( 5, *)DGAIN 

C set  header  block  for  ils  compatibility 
READd  ) (ICOMMONd  ) , 1 = 1,128) 

WRITE ( 2 ,REC=IREC) ( I COMMON ( I ) , 1=1,128) 

BLOCK=0 

250  READd , END=900 ) (ITRY( I ) , 1=1,256) 

BLOCK=l+BLOCK 
IF(BLOCK  .EQ.  1 ) THEN 
DO  141  J=l,256 
TTRY( 256+J)=REAL( ITRY( J) ) 

141  CONTINUE 
GO  TO  250 

ELSE 

IF(BLOCK  .EQ.  2 ) THEN 
DO  142  J=l,256 
TTRY( 512+J )=REAL( ITRY( J ) ) 

142  CONTINUE 
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C compute  SPDWD  kernel  for  first  block=0 
DO  143  J=257,512 
RIP=0 . 

DO  144  M=257-J,SMV 

RIP=TTRY(M+J)*TTRY(M+768-J)/(RSMV*DGAIN)+RIP 
144  CONTINUE 

IOUT( J-256 )=NINT( RIP ) 

143  CONTINUE 

IREC=IREC+1 

WRITE( 2,REC=IREC) ( IOUT( I ) , 1=1,256) 

GO  TO  250 
ELSE 

C generate  DWD  inner  product 
DO  200  J=l,256 
TTRY( J )=TTRY( 256+J ) 

TTRY( 256+J)=TTRY( 512+J) 

TTRY( 512+J )=REAL( ITRY( J) ) 

200  CONTINUE 

DO  211  J=257,512 
RIP=0 . 

DO  212  M=-SMV,SMV 

RIP=TTRY(M+J) *TTRY(M+768-J)/( RSMV*2000 . )+RIP 
212  CONTINUE 

IOUT( J-256 )=NINT( RIP) 

211  CONTINUE 

IREC=IREC+1 

WRITE ( 2 ,REC=IREC) ( IOUT( I ) , 1=1,256) 

GO  TO  250 
END  IF 
END  IF 

C compute  for  last  block 
900  DO  341  J=l,256 

TTRY( J)=TTRY( 256+J) 

TTRY( 256+J )=TTRY( 512+J) 

341  CONTINUE 

DO  342  J=257,512 
RIP=0 . 

DO  343  M=-SMV,512-J 

RIP=TTRY(M+J ) *TTRY(M+768-J)/(RSMV*DGAIN)+RIP 
343  CONTINUE 

IOUT( J-256 )=NINT( RIP) 

342  CONTINUE 
IREC=IREC+1 

WRITE( 2 ,REC=IREC) (ITRY(I)  , 1 = 1,256) 

CLOSE( 1 ) 

CLOSE( 2 ) 

STOP 

END 
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C CHRP2 

C This  program  generates  dual  chirp  data,  then  makes  it  ILS 
compatible 

INTEGER*4  ICOMMON(128) 

INTEGER*2  ITRY(256) 

REAL  F0,F1,TPI,RNU 
INTEGER  M 

CHARACTER*8  NAMEI,NAMEO 
PRINT  *,  'enter  FO' 

READ  *,  FO 

PRINT  *,  'enter  Fl' 

READ  *,  Fl 
WRITE( 6 ,100 ) 

100  FORMAT(lX,  'enter  the  sampled  data  file') 

READ( 5 , 110 )NAMEI 
110  FORMAT (A8) 

WRITE(6,120) 

120  FORMAT(lX,  'enter  the  chirp  file  name') 

READ( 5, 130 )NAMEO 
130  FORMAT (A8) 

OPEN( UNIT=1 ,NAME=NAMEI , STATUS= 'OLD' , FORM= ' UNFORMATTED' 

1 READONLY) 

OPEN ( UNI T=2 , NAME=NAMEO , STATUS= ' NEW' , FORM= ' UNFORMATTED ' , 

1 ACCESS='DIRECT' ,RECL=128 ) 

C initialize  parameters 

PRINT  *, 'enter  data  pnts  per  period' 

READ  * , RNU 
TPI=8 . *ATAN( 1 . ) 

IREC=1 

M=-256 

C set  header  block  for  ils  compatibility 
READ( 1 ) ( ICOMMON( I ) , 1=1,128) 

WRITE( 2 ,REC=IREC) ( I COMMON ( I ) , 1=1,128) 

250  READ( 1 ,END=900 ) ( ITRY( I ) , 1=1,256) 

C Generate  chirp  data 
M=M+256 

DO  200  J=l,256 

ITRY( J )=1000 . *COS ( TPI*F0* ( ( M+J ) /RNU) **2/2 . ) 

1 +1000.*COS(TPI*F1*( (M+J)/RNU)**2/2. ) 

200  CONTINUE 

IREC=IREC+1 

WRITE( 2 ,REC=IREC) (ITRY(I) , 1=1,256) 

GO  TO  250 
900  CONTINUE 

CLOSE( 1 ) 

CLOSE( 2 ) 

STOP 

END 
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C ECH02D 

generates  data  for  two  time  displaced  FM 

^ then  makes  it  ILS  compatible 

INTEGER*4  ICOMMON(128)  pdcioie 

INTEGER*2  ITRY( 256 ) , ITRY2 ( 256 ) 

REAL  F0,TPI ,RNU,GAIN,GAINE 
INTEGER  M,DISP 
CHARACTER*8  NAMEI, NAMED 
PRINT  *,  'enter  FO' 

READ  *,  FO 

PRINT  *,  'enter  Fl ' 

READ  *,  Fl 

WRITE(6,100) 

sampled  data  file') 

READ(5,110)NAMEI 
FORMAT (A8 ) 

WRITE(6,120) 

FORMAT(lX,  'enter  the  chirp  file  name') 
READ(5,130)NAMEO  ^ 

FORMAT (A8 ) 


100 

110 

120 

130 


0PEN(UNIT=1^NAME=NAMEI,STATUS='0LD' ,FORM='UNFORMATTED' , 

OPEN ( UNI T=2 , NAME=NAMEO , STATUS= ' NEW ' , FORM= ' UNFORMATTED ' 

1 ACCESS='DIRECT' ,RECL=128)  wKnAiXbu  , 

C initialize  parameters 

1536  )'  *'  'enter  data  points  per  period  (eg 

READ  *,RNU 

PRINT  *,  'enter  time  displacement  value' 

READ  *,  DISP 

PRINT  *,  'enter  the  signal  amplitude' 

READ  *,  GAIN 

PRINT  *,  'enter  the  echo  ampl ' 

READ  *,  GAINE 
TPI=8 . *ATAN( 1 . ) 

IREC=1 

M=-256 

C set  header  block  for  ils  compatibility 
READ(l)(ICOMMON(I),  1=1,128) 

WRITE(2,REC=IREC) (ICOMMON(I) , 1=1,128) 

250  READ( 1 ,END=900 ) ( ITRY( I ) , 1=1,256) 

C Generate  chirp  data 
M=M+256 
M2=M-DISP 
DO  200  J=l,256 

ITRY( J ) =GAIN*COS (TPI*F0*( (M+J) /RNU ) **2/2 . ) 

200  ^JTRY2U)=gaINE*COS(TPI*F1*( (M2+J)/RNU)**2/2. ) 

DO  300  J=l,256 


ITRY( J)  = ITRY( J)  + ITRY2(  J) 

300  CONTINUE 

IREC=IREC+1 

WRITE( 2,REC=IREC) ( ITRY( I ) , 1=1,256) 
GO  TO  250 
900  CONTINUE 

CLOSE( 1 ) 

CLOSE( 2 ) 

STOP 

END 
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C ECHO 

C This  program  generates  data  for  two  time  displaced  FM 
signals , 

C then  makes  it  ILS  compatible 
INTEGER*4  ICOMMON(128) 

INTEGER*2  ITRY( 256 ) , ITRY2 ( 256 ) 

REAL  F0,TPI,RNU 
INTEGER  M,DISP 
CHARACTER*8  NAMEI,NAMEO 
PRINT  *,  'enter  FO' 

READ  *,  FO 
WRITE(6,100) 

100  FORMAT(lX,  'enter  the  sampled  data  file') 

READ( 5 , 110 )NAMEI 
110  FORMAT (A8) 

WRITE(6,120) 

120  FORMAT(lX,  'enter  the  chirp  file  name') 

READ( 5,130)NAMEO 
130  FORMAT (A8) 

OPEN( UNIT=1 , NAME=NAMEI , STATUS= ' OLD' , FORM= ' UNFORMATTED' , 
1 READONLY) 

OPEN( UNIT=2 , NAME=NAMEO , STATUS= ' NEW' , FORM= ' UNFORMATTED' , 
1 ACCESS='DIRECT' ,RECL=128 ) 

C. initialize  parameters 

PRINT  *,  'enter  data  points  per  period  (eg 

1536.)' 

READ  *,RNU 

PRINT  *,  'enter  time  displacement  value' 

READ  *,  DISP 
TPI=8 . *ATAN( 1 . ) 

IREC=1 

M=-256 

C set  header  block  for  ils  compatibility 
READ( 1 ) ( ICOMMON( I ) , 1=1,128) 

WRITE( 2,REC=IREC) ( ICOMMON( I ) , 1=1,128) 

250  READ(1,END=900) (ITRY(I) , 1=1,256) 

C Generate  chirp  data 
M=M+256 
M2=M-DISP 
DO  200  J=l,256 

ITRY( J ) =1000  . *COS ( TPI*F0* ( ( M+ J ) /RNU ) * * 2/2 . ) 
ITRY2 ( J )=1000 , *COS ( TPI*F0* ( ( M2+ J ) /RNU ) * *2/2 . ) 
200  CONTINUE 

DO  300  J=l,256 
ITRY( J)=ITRY( J)+ITRY2( J) 

300  CONTINUE 

IREC=IREC+1 

WRITE( 2 ,REC=IREC ) ( ITRY( I ) , 1=1,256) 

GO  TO  250 
900  CONTINUE 


204 


CLOSE(l) 
CLOSE( 2 ) 
STOP 
END 
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C SPDWDW2 

C This  program  takes  ILS  compatible  sampled  data, 

C windows  it,  generates  the  SPDWD 

C inner  product,  then  outputs  the  IP  as  an  ILS  compatible 
data  file; 

C in  this  form,  when  passed  through  the  ILS  program  SDI  the 
output  is 

C the  SPDWD  of  the  original  data. 

INTEGER*4  ICOMMON(128) 

INTEGER*2  ITRY( 256 ) , IOUT( 256 ) , SMV 

REAL  TTRY(768) , RIP , FLTR( 256 ) ,WNDW(256) ,WNDWS(256) 

CHARACTER*8  NAMEI, NAMED 

WRITE( 6 ,100 ) 

100  FORMATdX,  'enter  the  sampled  data  file') 

READ( 5 , 110 )NAMEI 
110  FORMAT (A8) 

WRITE( 6 , 120 ) 

120  FORMATdX,  'enter  the  SPDWD  IP  file  name') 

READ( 5 , 130 ) NAMED 
130  FORMAT (A8) 

OPEN(UNIT=l ,NAME=NAMEI , STATUS= ' OLD ' , FORM= ' UNFORMATTED ' , 

1 READONLY) 

OPEN ( UNIT=2 , NAME=NAMEO , STATUS= ' NEW ' , FORM= ' UNFORMATTED ' , 

1 ACCESS='DIRECT' ,RECL=128 ) 

C initialize  parameters 
IREC=1 
WRITE( 6 , 11 ) 

11  FORMAT(lX,  'enter  the  smoothing  parameter  value') 

READ(5,*)SMV 

RSMV=REAL( SMV) 

WRITE( 6 , 12 ) 

12  FORMAT(lX,  'enter  the  necessary  magscale  for  ILS 
compatibility' ) 

READ( 5,*)DGAIN 

C generate  window  and  filter  functions 
READ( 3,*) (WNDW( J) , J=l,257) 

DO  33  J=l,257 

WNDWS ( J ) =WNDW ( J ) *WNDW ( 2 5 7 - J ) 

33  CONTINUE 

PRINT  *,  'enter  the  Gaussian  attenuation  parameter' 
READ  *,  ALPHA 
DO  34  J=1,SMV 
R=-l . *ALPHA*J**2 
FLTR( J )=EXP(R) 

34  CONTINUE 

C set  header  block  for  ils  compatibility 
READd  ) (ICOMMONd  ) , 1 = 1,128) 

WRITE( 2,REC=IREC) ( ICOMMONd ) , 1=1,128) 

BLOCK=0 

READd  ,END=900  ) (I  TRY  (I  ) , 1 = 1,2  56) 


250 


BL0CK=1+BL0CK 
IF(BLOCK  .EQ.  1 ) THEN 
DO  141  J=l,256 
TTRY(256+J)=REAL(ITRY( J) ) 

141  CONTINUE 
GO  TO  250 

ELSE 

IF(BLOCK  .EQ.  2)THEN 
DO  142  J=l,256 
TTRY( 512+J)=REAL( ITRY( J) ) 

142  CONTINUE 

C compute  SPDWD  kernel  for  first  block=0 
DO  143  J=257,512 
RIP=0 . 

DO  144  M=257-J,SMV 
IF(M  .LT.  0)THEN 
MP=-1*M 
ELSE 
MP=M 
END  IF 

IF(MP  .EQ.  0)THEN 
RIP=RIP 
ELSE 

RIP=TTRY ( M+J ) *TTRY( M+768-J ) *FLTR( MP ) /( RSMV*DGAIN ) +RI 
END  IF 

144  CONTINUE 

OUT=RIP*WNDWS ( J-256 ) 

IOUT( J-256 )=NINT( OUT) 

143  CONTINUE 
IREC=IREC+1 

WRITE( 2,REC=IREC) (lOUT(I) , 1=1,256) 

GO  TO  250 
ELSE 

C generate  DWD  inner  product 
DO  200  J=l,256 
TTRY( J)=TTRY( 256+J) 

TTRY( 256+J)=TTRY( 512+J) 

TTRY( 512+J )=REAL( ITRY( J ) ) 

200  CONTINUE 

DO  211  J=257,512 
RIP=0. 

DO  212  M=-SMV,SMV 
IF(M  .LT.  0)THEN 
MP=-1*M 
ELSE 
MP=M 
END  IF 
IF(MP  .EQ. 

RIP=RIP 
ELSE 


0 )THEN 
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RIP=TTRY(M+J ) *TTRY( M+768-J ) *FLTR( MP )/( RSMV*2000 . )+RIP 


END  IF 

212 

CONTINUE 

OUT=RIP*WNDWS( J-256 ) 
IOUT( J-256 )=NINT( OUT) 

211 

CONTINUE 

IREC=IREC+1 

WRITE( 2 ,REC=IREC) ( IOUT( I ) 
GO  TO  250 
END  IF 
END  IF 

C compute 

for  last  block 

900 

DO  341  J=l,256 
TTRY ( J ) =TTRY ( 2 56  + J ) 
TTRY( 256+J)=TTRY( 512+J) 

341 

CONTINUE 

DO  342  J=257,512 
RIP=0 . 

DO  343  M=-SMV,512-J 
IF(M  .LT.  0)THEN 
MP=-1*M 
ELSE 
MP=M 
END  IF 

IF(MP  .EQ.  0)THEN 
RIP=RIP 
ELSE 

RIP=TTRY( M+J ) *TTRY( M+768-J ) *FLTR( MP )/( RSMV*DGAIN ) +RI P 
END  IF 

343  CONTINUE 

OUT=RIP*WNDWS( J-256 ) 

IOUT( J-256 )=NINT(OUT) 

342  CONTINUE 

IREC=IREC+1 

WRITE( 2,REC=IREC) ( ITRY( I ) , 1=1,256) 

CLOSE( 1 ) 

CLOSE( 2 ) 

STOP 

END 
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C SPLUSN 

C This  program  adds  noise  to  data,  then  makes  it  ILS 
compatible 

INTEGER*4  ICOMMON( 128 ) ,M,NRIP 
INTEGER*2  ITRY(256) 

REAL  FO,TPI,RNU,RIP,THR 
CHARACTER*8  NAMEI,NAMEO 
WRITE(6,100) 

100  FORMATdX,  'enter  the  sampled  data  file') 

READ ( 5,110 )NAMEI 
110  FORMAT (A8) 

WRITE(6,120) 

120  FORMAT(lX,  'enter  the  signal  plus  noise  file  name') 

READ(5,130)NAMEO 
130  FORMAT (A8) 

OPEN ( UNI T=1 ,NAME=NAMEI , STATUS= ' OLD' , FORM= ' UNFORMATTED ' 

1 READONLY) 

OPEN ( UNIT=2 ,NAME=NAMEO, STATUS= ' NEW' , FORM= ' UNFORMATTED ' , 

1 ACCESS='DIRECT' ,RECL=128 ) 

C initialize  parameters 

PRINT  *,  'enter  seed' 

READ  *,M 

PRINT  *,  'enter  noise  threshold' 

READ  *,THR 
IREC=1 

C set  header  block  for  ils  compatibility 
READ( 1 ) ( ICOMMON( I ) , 1=1,128) 

WRITE( 2 ,REC=IREC) ( ICOMMON( I ) , 1=1,128) 

250  READ( 1 , END=900 ) ( ITRY( I ) , 1=1,256) 

C Generate  noise  data 

DO  200  J=l,256 
RIP=RAN(M) 

RIP=THR*RIP 

NRIP=NINT(RIP) 

ITRY( J)=ITRY( J)+NRIP 
200  CONTINUE 

IREC=IREC+1 

WRITE( 2 ,REC=IREC) ( ITRY( I ) , 1=1,256) 

GO  TO  250 
900  CONTINUE 

CLOSE( 1 ) 

CLOSE( 2 ) 

STOP 

END 
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C 

C 


100 


SPNGAU 

This  program  adds  Gaussian  noise  to  ILS  data 
INTEGER*4  ICOMMON ( 1 2 8 ) , M , NRI P 
INTEGER*2  ITRY( 256 ) , IGAU( 256 ) 

REAL  F0,TPI,RNU,RIP,STD,RMEAN 
CHARACTER*8  NAMEI , NAMEO , NAMEN 
WRITE(6,100) 

FORMATdx,  'enter  the  sampled  data  file') 
READ(5,110)NAMEI  ’ 


110  FORMAT(A8) 

WRITE(6,120) 

120  FORMATdx,  'enter 

READ( 5,130 ) NAMEO 
130  FORMAT(A8) 

WRITE(6,140) 

140  FORMATdx,  'enter 

READ( 5 , 150 ) NAMEN 
150  FORMAT (A8) 


the  signal  plus  noise  file  name') 


the  noise  file  name') 


OPEN(UNIT=dNAME=NAMEI,STATUS='OLD' ,FORM='UNFORMATTED' , 

OPEN ( UNI T=2 , NAME=NAMEO , STATUS= ' NEW ' , FORM= ' UNFORMATTED ' 

1 ACCESS='DIRECT' ,RECL=128)  , 

OPEN ( UNI T= 3 , NAME=NAMEN , STATUS= ' NEW ' , FORM= ' UNFORMATTED ' 

1 ACCESS='DIRECT' ,RECL=128) 

C initialize  parameters 

PRINT  *,  'enter  seed' 

READ  *,M 

PRINT  *,  'enter  seed  coefficient' 

READ  *,  GAIN 

PRINT  *,  'enter  the  large  number' 

READ  * , LNU 
HLNU=LNU/2 . 

PRINT  *,  'enter  the  standard  deviation' 

READ  *,  STD 

PRINT  *,  'enter  the  mean' 

READ  * , RMEAN 
IREC=1 

C set  header  blocic  for  ils  compatibility 
READd  ) (iCOMMONd  ) , 1 = 1,128) 

WRITE(2,REC=IREC) dCOMMON(I) , 1=1,128) 

WRITE( 3, REC=IREC) (ICOMMONd ) , 1=1,128) 

250  READd,  END=900  ) (ITRYd)  , 1 = 1,256  ) 

C Generate  Gaussian  noise  via  law  of  large  numbers 
DO  199  1=1,256 
GA=0  . 

DO  200  J=1,LNU 
RIP=RAN(M) 

RM=RIP*GAIN 
M=NINT( RM) 
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200 

199 

210 


900 


GA=RIP+GA 

CONTINUE 

RIP2= ( GA-HLNU ) *STD+RMEAN 
IGAU( I )=NINT(RIP2 ) 

CONTINUE 
DO  210  J=l,256 
ITRY( J)=ITRY( J)+IGAU( J) 

CONTINUE 

IREC=IREC+1 

WRITE( 2 ,REC=IREC) ( ITRY( I ) , 1=1,256) 
WRITE( 3,REC=IREC) ( IGAU( I ) , 1=1,256 ) 
GO  TO  250 
CONTINUE 
CLOSE( 1 ) 

CLOSE( 2 ) 

CLOSE( 3 ) 

STOP 

END 
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C RFILT 

C This  program  takes  sampled  data  and  filters  it  via 

C a recursive  filter 

C 


INTEGER*4  ICOMMON(128) 

INTEGER*2  ITRY(256) 

DIMENSION  ZLOW(9) ,ZBAND(9) , AL ( 9 ) , PL ( 9 ) , AB ( 9 ) , PB ( 9 
DIMENSION  CC(IO) 

CHARACTER*8  NAMEI,NAMEO 
CHARACTER* 8 NAMEA,NAMEP 
WRITE(6,100) 

FORMATdX,  'Enter  the  sampled  data  file  to  be 
processed ' ) 

READ( 5 , 110 )NAMEI 
FORMAT (A8 ) 

WRITE(6,120) 

FORMATdX,  'Enter  the  processed  data  file  name') 
READ( 5 , 130 )NAMEO 
FORMAT (A8 ) 

WRITE(6,140) 

FORMAT(lX,  'Enter  the  denom  coeffs') 

READ( 5 , 150 )NAMEA 
FORMAT (A8 ) 

WRITE(6,160) 

FORMAT(lX,  'Enter  the  num  coeffs') 

READ( 5,170)NAMEP 
FORMAT (A8 ) 

OPEN 

( UNIT=1 , NAME=NAMEI , STATUS= ' OLD ' , FORM= ' UNFORMATTED ' 

1 READONLY) 

OPEN 


100 

pro 

110 

120 

130 

140 

150 

160 

170 


( UNIT-2 , NAME=NAMEO , STATUS= ' NEW' , FORM= ' UNFORMATTED ' 

1 ACCESS='DIRECT' ,RECL=128 ) 

OPEN  (UNIT=3,NAME=NAMEA,STATUS='OLD' , READONLY) 
^ OPEN  ( UNIT=4,NAME=NAMEP,STATUS=' OLD ', READONLY) 

C initialize  params  and  read  filter  coeffs 


ICOSI=0 

IREC=1 

READ(3,*)  (AB(I),  1=1,9) 

READ(4,*)  (PB(I),  1=1,9) 

C 

C set  header  block  for  new  data  file 

C 


READ(l)  dCOMMON(I),  1=1,128) 

WRITE( 2,REC=IREC)  (ICOMMON(I),  1=1,128) 
C 

550  READd , END=900 ) (ITRY(I),  1=1,256) 

C 

C 

C Filter  loop 


o o o 
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C 

C 


900 


DO  L=l,256 


TIM=0 . 

DO  IBP=2,9 

ZBAND( 11-IBP)=ZBAND( 10-IBP) 
TIM=TIM-AB( 11-IBP) *ZBAND( 11-IBP) 
END  DO 

ZBAND( 1 )=TIM+FLOAT( ITRY( L ) ) 

YI  = 0. 

DO  IBPF=1,9 

YI=YI+ZBAND( IBPF)*PB( IBPF) 

END  DO 

gain 

ITRY(L)=NINT(  YI ) 

END  DO 
IREC=IREC+1 

WRITE(2,REC=IREC)  (ITRY(I),  1=1, 
GOTO  550 
CONTINUE 
CLOSE( 1 ) 

CLOSE( 2 ) 

CLOSE( 3 ) 

CLOSE( 4 ) 

STOP 

END 


256) 


C DEMZOOM 

C this  program  modulates  ILS  compatible  data 
INTEGER*4  ICOMMON(128) 

INTEGER*2  ITRY(256) 

REAL  FC,TPI ,RNU,RN 
INTEGER  M 

CHARACTER*8  NAMEI, NAMED 
PRINT  *,  'enter  FC' 

READ  * , FC 

PRINT  *,  'enter  period' 

READ  * , RNU 
WRITE( 6 ,100 ) 


100 

FORMATdX,  'enter 
READ( 5 , 110 )NAMEI 

the 

sampled  data  file') 

110 

FORMAT (A8 ) 
WRITE( 6 , 120 ) 

120 

FORMATdX,  'enter 
READ( 5 , 130 )NAMEO 

the 

processed  file  name') 

130 

FORMAT (A8 ) 

OPEN{ UNIT=1 , NAME=NAMEI , STATUS= ' OLD' , FORM= ' UNFORMATTED 
1 READONLY) 

OPEN( UNIT=2 ,NAME=NAMEO, STATUS= 'NEW' , FORM= ' UNFORMATTED 
1 ACCESS='DIRECT' ,RECL=128 ) 

C initialize  parameters 
TPI=8 . *ATAN( 1 . ) 

IREC=1 

M=-256 

C set  header  block  for  ils  compatibility 
READ( 1 ) ( ICOMMON( I ) , 1=1,128) 
WRITE(2,REC=IREC) ( I COMMON ( I) , 1=1,128) 

250  READ( 1 , END=900 ) ( ITRY( I ) , 1=1,256) 

C heterodyning  process 
M=M+256 

DO  200  J=l,256 

ITRY( J )=ITRY( J ) *COS ( TPI*FC* ( M+J ) /RNU) 

200  CONTINUE 

IREC=IREC+1 

WRITE(2,REC=IREC) (ITRY(I) , 1=1,256) 

GO  TO  250 
900  CONTINUE 

CLOSE( 1 ) 

CLOSE( 2 ) 

STOP 

END 
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